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Chapter 5 


Many Body Perturbation Theory 


Density Functional Theroy (DFT) calculations are sufficiently accurate for many applications. 
However its description of electronic correlation is just approximated locally from the uniform 
electron gas. It lacks, for instance, indirect effects of the electron-electron interaction, such 
as Van-der-Walls force where electrostatic repulsion deforms the distribution of electrons and 
these polarized electron distributions in turn can attract each other. This effect can not be 
captured by DFT or HF as they only take the instantaneous electrostatic interaction into 
consideration but the formation of the polarized electron distributions requires time. The 
missing descriptions of such dynamical effects is referred to as dynamical correlation (Shavitt 
and Bartlett 2009). 

Static correlation occurs when the DFT or HF approximation is not appropriate for the 
chemical environment per se. A good example is the the dissociation of a hydrogen ion where 
the protons are already at a great distance. The single electron should be in a superposition 
of two states, quite localized at each respective proton. DFT rather places half an electron 
on each proton and yields just one orbital spanning both protons. 

Perturbation theory expands properties of the exact solution of the Schrodinger equation 
in terms of the orbitals and orbital energies of the corresponding Hartree-Fock approximation. 
One can also start from DFT orbitals or from other reference which is feasible to solve. It 
is mostly the dynamic correlation that is captured by the perturbation expansion but it can 
also include some static correlation when going to sufficiently many terms. In cases where the 
static correlation is large multi reference perturbation theory may be required. Not unlike a 
Taylor expansion, a perturbation expansion is not guaranteed to converge or may converge 
slowly with the number of terms included. The convergence also strongly depends on the 
quality of the reference state. 

We will use time dependent many body perturbation theory following the original deriva¬ 
tion of (Goldstone 1957) as it is independent of the reference system (DFT or HF) and it 
is extensive. Extensivity means that for two systems A and B the energy of the combined 
system equals the sum of the energy of the individual constituents 

E{AB) = E{A) + E[B), 

assuming an identical chemical environment. The time dependent formulation of perturbation 
theory also lends itself naturally to Goldstone diagrams to visualize terms occurring in the 
perturbation expansion. 

Complementary treatments can be found in (Kutzelnigg 2009; Lancaster and Blundell 
2014; Shavitt and Bartlett 2009; Fetter and Walecka 2003; Coleman 2015). 
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CHAPTER 5. MANY BODY PERTURBATION THEORY 


5.1 Time dependent perturbation theory 

In perturbation theory the exact Hamiltonian is separated into the unperturbed part Hq, 
which can be solved, and the perturbation Hi, which contains the full Coulomb electron- 
electron interaction: 

H = T -\- Vne + YeS p Vee ~ Yq . 

' -V-' '-" 

i/o ill 

I4ff is the effective interaction employed by the reference, e.g. the Hartree-Fock approxima¬ 
tion. Note that this effective interaction included in Hq must be subtracted again by the 
perturbation to arrive at results of the full Hamiltonian. 

Let V'p be the spin-orbitals of the unperturbed Hamiltonian of the HF or DFT reference 
and let |<h) denote the Slater determinant of the ground state, where the lowest N states are 
occupied by the N electrons present in the system. These states are called unexcited states 
while the states that are unoccupied in |$) are called excited states. We will use the letters 
i,j,k,... to label unexcited states, a,b,c,... to label excited states and p,q,r,... to label 
general states. We can write |<I>) in second quantization as the result of applying the electron 
creation operator ct for all unexcited states i on the vacuum state, | ), without any electrons: 

1^)=n^ii) 

i 

Note that each application of c| changes the number of particles and thus the dimensionality 
of the Hilbert space. The states of second quantization are elements of the union of all Hilbert 
spaces of zero particles, one particle, two particle and so forth, which is called Fock-space. The 
beauty of second quantization is that it hides all the tedious footwork of anti-symmetrization 
in the algebra of the creation and annihilation operators, which is completely given by the 
anti-commutator relations: 

Yl, c\} = 0 {Cp, Cq} = 0 {cj,, Cq] = 5pq, (5.1) 

where {A,B} = AB + BA. One immediate consequence of the these relations is the Pauli 
exclusion principle disallowing two fermions in the same state: 

I ) ^ 

Note that the vacuum, | ), is a state while the number 0 is not. 

The Fock-space used for the creation and annihilation operators is spanned by the Slater 
determinants of the eigenfunctions ijjp of the unperturbed Hamiltonian Hq. Hq is therefore 
diagonal, counting the eigenenergy £p for each occupied state p, irrespective of whether it is 
an excited or unexcited state: 

Hq = ^ ^ ^pCpCp. (5-2) 

p 


5.1.1 Particle/hole picture 

Let us now introduce the particle/hole picture where we want to consider the non-interacting 
ground state |<I>) as the new vacuum state instead of the true vacuum, | ), without any 
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electrons. In this picture we care about the difference to the non-interacting ground state |<I>) 
and we only count excited states that are now occupied called particles, and unexcited states 
that are no longer occupied called holes. While cl indeed creates a particle in an excited 
state, a hole in an unexcited state has to be created by annihilating a formerly occupied state 
by Cj. For unexcited states below the Fermi energy the meaning of creation and annihilation 
has to be reversed. We define the creation and annihilation operators for particles (a^,a) and 
for holes (i^, i): 

= cl a = Ca for 2 ^ 

i^ = Ci 'I = cl for Ei < ep 

The anti-commutator relations for these operators follow directly from (5.1) and the only 
non-vanishing relations are 

{i\j} = dij {a\b} = dab- (5.4) 

Inserting the particle/hole operators into (5.2) splits the sum over states p into a sum over 
holes i and a sum over particles a, giving 


hq = Eiii^ + 






SaOj 


where we have used the anti-commutator relation to put the operators into normal order, 
such that all annihilation operators appear on the right. 

Next, we need to translate the perturbation Hi into the particle/hole formalism. We 
start with its second quantized representation using the electron creation and annihilation 
operators Cp and c^: 


pqrs 




^C^clcq, 


where Vir and Vq are the matrix elements of the first quantized operators I4e and I4fr) given 
by 


yir = {vq\yee\sr) = / / dx dx'V’* (x) V’* (x) ^^ V’r (x) V's (x) 


< = (p\yes\q) = / dxV’p(x) 


with f dx = I br. The factor 4 in (5.6) accounts for double counting when not restricting 
the sum over p, q, r and s to distinct elements of V^r. 

Unfortunately, it is not as straight forward to translate Hi into the particle/hole picture 
as it was for Hq since there are now up to 4 different indices which can be either holes or 
particles. If, for example, p, q, r are particle indices a, c, b and s is a hole index i, Vee creates 
the particles a, c and the hole i while it destroys the particle b, as shown in Figure 5.1. The 
electron creation and annihilation operators and c always occur in pairs thus keeping the 
number of electrons constant. However, particle and hole creation and annihilation operators 
are in general not normal ordered and occur in any constellation. The number of particles 
and holes is therefore not necessarily constant. Figure 5.1 also shows a particle/hole pair 
created by the action of Ugg indicated by a shaded circle. 
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Figure 5.1: Examples of terms occurring in Hi and their respective diagrammatic represen¬ 
tation 

5.1.2 Interaction pictnre 

So far, all operators were given in the Schrodinger picture, where all time evolution takes 
place in the states and the operators are time independent. The Heisenberg picture, on 
the other hand, keeps the states time independent and all time evolution is put into the 
operators. In order to do time dependent perturbation theory, it is convenient to use the 
interaction picture, which is a hybrid of the Schrodinger and the Heisenberg picture. In the 
interaction picture the operators evolve according to a time evolution solely based on the 
unperturbed part of the Hamiltonian Hq, while the states evolve only due to the action of 
the perturbation Hi. 

Given a state in the Schrodinger picture |'k(t)) we define the corresponding state in the 
interaction picture |T 7 (t)) by “undoing” the time evolution which originates from Hq: 

= (5.9) 

To find the equation of motion for states in the interaction picture we derive with respect to 
time, giving 

ik |5.,(()) = |<P(()) + |4.(t)> 

|5.{()) + + Hi) |4'(t)> 

= P^°^Hie-^^°^\'Ifi{t)). (5.10) 

'-v^ 


So the time evolution of |T/(t)) is determined by the perturbation Hii{t) only. The time 
dependence of Hii{t), on the other hand, solely depends on the unperturbed part of the 
Hamiltonian Hq. We are now interested in the time evolution operator Ui{t,to) that evolves 
a state in the interaction picture from the time to to the time t: 

Ui{t,to)\^i{to)) = \^i{t)). (5.11) 

Operating with id/dt on both sides and using (5.10) gives 

i-^Ui{t,to) |T/(to)) = Hii{t)Ui{t,to) |'k/(to)) , 

which must holds for all |T/(to)) thus leading to the equation of motion for the time evolution 
operator in the interaction picture 

Q .... 

i—Ui{t,to) = Hij{t)Ui{t,to). 


(5.12) 
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Integrating both sides of the above equation with respect to the first argument from to to t 
yields 

Ui{t, to) = 1 - i / dt' Hu{t')Ui{t', to), (5.13) 

Jto 

where we have used Ui{to,to) = 1, which follows from (5.11). This integral equation can be 
solved by iteratively applying the above equation to each occurrence of Uj{t',to),Ui{t”,to) 
and so forth, giving 

rt p'b 

Ui{t,to) = 1-i/ dt'Hii{t') + \^ dt' dt" Hu{t')Hu{t")--i^... 

J tQ J tQ J tQ 

OO « 

= / dH...dtnHii{ti)...Hii{tn). (5.14) 

This means that the time evolution operator in the interaction picture can be constructed 
by applying the perturbation Hu any number of times and at all possible times between to 
and t in a time ordered manner. The number of applications of Hu is called the order in the 
perturbation expansion. The perturbation expansion for the time evolution operator can be 
truncated at a finite order. However, the truncated expansion is not guaranteed to converge 
at any finite order. In metals, for example, no truncation is convergent beyond first order. 
This is a considerable drawback for finite order perturbation methods, such as second order 
Mpller-Plesset perturbation theory (MP2). 

5.2 The Gell-Mann—Low theorem 

In the previous section we have seen how to evolve any given state from to to t. We are, 
however, interested in eigenstates of the full Hamiltonian, most importantly in its ground 
state jT). All we have are eigenstates of the non-interacting Hamiltonian Ho which we can 
evolve from. In general, such an evolved state Ui{t,to) |‘I'(to)) will be neither an eigenstate 
of the non-interacting Hamiltonian Ho nor of the full Hamiltonian H. We can, however, 
introduce a time dependent perturbation 

H{t) = Ho + e'>^Hi (5.15) 

for t < 0, slowly turning on the electron-electron interaction for small rj > 0. The system 
starts with the unperturbed Hamiltonian Hq at t = — oo and at t = 0 the interaction is fully 
turned on, giving H. According to the adiabatic theorem, the system will stay in an eigenstate 
of H{t) at all times when starting from an eigenstate of Ho at t = —oo if the transition is 
sufficiently slow, which holds for 

rj <C AE. 

AE is the energy change of the ground state when turning on the interaction. Note that 
starting in the ground state of Ho at t = —oo does not guaranteed that the evolved eigenstate 
at t = 0 is indeed the ground state of the full Hamiltonian H. Level crossings may occur. We 
will, however, assume that they do not occur for the ground state and that the ground state 
|<h) of Ho evolves adiabatically into the ground state I'k) of H: 


|T) = ^^(0,-oo)|<h) 


( 5 . 16 ) 
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where we drop the explicit denotation of the interaction picture in favor of denoting the 
dependence on the parameter rj of the transition speed, writing from now on 

oo „ 

Ur,{t,to) = Y.{-ir / dh...dtnHl{tl)...Hi{tn) (5.17) 

t>tl>...>t„>tQ 


and 

Fi(t) = (5.18) 

Although Urj and Hi depend on r] we hope that the results are, in the end, independent of 
T] if it is chosen sufficiently small. The quantity we are now most interested in is the energy 
difference AE between the ground state energy of the non-interacting Hamiltonian Hq and the 
ground state energy of the fully interacting Hamiltonian H = Hq + Hi. The fully interacting 
ground state |T) satisfies the Schrodinger equation 

{Hii + Hi)\^>) = {Eo + AE)\^>). 


Multiplying both sides with (<1>| from the left and using (5.16) gives 


The last equation now relates the ground state energy of the fully interacting system to 
vacuum expectation values (VEV) of the non-interacting system in the particle/hole picture, 
where |<1>) is the vacuum state without any particle or hole excitations. In principle, these 
expectation values can be evaluated despite the infinite sums hidden in the time evolution 
operators tin- 


(^|i7i|T) 

(d>|T) 

-oo)|$) ^ (d>|ffi(0)17^(0,-oo)|^) 
> 0 ) 1 $) {mr,{0,-oo)m 


5.3 Wick’s theorem 


Equation (5.19) still has two practical drawbacks. First, there is yet no systematic recipe 
given how to approximate the time evolution operator Uj), consisting of an infinite sum, and 
second, one still needs to show that AE does not depend on the choice of rj if it is chosen 
sufficiently small. We start with the diagrammatic representation of the terms occurring in 
(5.19), where we need to evaluate terms of the form 

A,B,... are arbitrary creation or annihilation operators. If these operators were normal 
ordered, such that all annihilation operators are on the right side of all creation operators, 
its vacuum expectation value (VEV) would simply be 0, given there is at least one creation 
or annihilation operator. A relation between the VEV of the operators in the order given 
and the VEV in normal order is therefore desired. Let N[AB ...] denote normal ordering of 
the operators where in the case of fermions the sign changes for each transposition. For two 
creation or annihilation operators, there is only one non-trivial case 

N\pq^ = —q^p, 
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in the other three cases the operators are already in normal order. Since the vacuum ex¬ 
pectation value of two operators in normal order is 0, we can rewrite a given VEV of two 
operators 


(^>|i5|^>) = - ($|7V[ii3]|^>) = {^\AB - 7V[i.B]|^>) = 

defining the contraction of two operators A and B by^ 


(5.20) 


AI3 := AB - N[Ab]. 

This may seem arbitrary but it turns out that the contraction of two operators is simply a 
number. For two creation or annihilation operators there are four possible contractions 


pq = 0 


pq 


t = gt} = pt^ = 0 p^q^ = 0. 


(5.21) 


Since the contraction of two operators is a number it is unaffected by normal ordering. This 
allows us to make the desired connection between the given order and the normal order of 
two operators AB: 

I—I I—I 

Ab = n[Ab] + Ab = n[Ab + Ab] 

This result can be generalized to more than two operators by the virtue of Wick’s theorem 
(Wick 1950; Peskin and Schroeder 1995): 


ABC ... = N 


ABC ... + all possible contractions of ABC ... 


(5.22) 


For a sequence of four operators this gives for instance 

I—I I 1 

AbCD = N[AbCD] -f N[AbCD] + N[AbCD] + N[ABCD] + N[ABCD] -f 


I-1 I —I I —II—I I I—I— I I I —I I 

n[Abcd] -f n[Abcd] + n[Abcd] + n[Abcd] + n[Abcd], ( 5 . 23 ) 

where we can now reorder the operators, changing the sign appropriately, to pull contracted 
operators out of the normal ordering operator. For example 


I-1 I—I 

N[AbCD] = -Ac X N[BD]. 

Since vacuum expectation values of normal ordered operators vanish, only the fully contracted 
terms survive. The VFV of the four operators in (5.23) thus evaluates to 


I—^ I—I I—^ I—I I—^ I—I 

{^\Abcd\<i>) = Abcd - Acbd + Adbc. ( 5 . 24 ) 


I—I 

^Note that contractions are usually defined by AB — T[AB] — A^[AB], where T[AB] denotes time ordering 
of the operators, such that the time increases from right to left. In our case, the terms occurring in (5.19) are 
already time ordered by the constraints of the integrals so we will drop the time ordering symbol. 
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5.3.1 Application to the perturbation 

We can now apply this algebra to a simplified case without an effective interaction where the 
perturbation is given by 


H,{t) 




1 

2 


'^yirClc\crCse . 

pqrs 


(5.25) 


Note that the electron creation and annihilation operators Cp and Cp have to be expressed in 
terms of particle and hole creation and annihilation operators and p, since we want to use 
the non-interacting ground state |<h) as the vacuum state. From (5.3) we get 



p^ for Sp > ep 
p otherwise, 


J p for £p > ep 
( p^ otherwise. 


(5.26) 


We start with evaluating the numerator of the right term in (5.19) in zeroth order of the 
expansion of the time evolution operator, where —oo) = 1 : 


(^l-^l(O)l^) = \Y1 VP^clc\crCs + VP^clc\crCs + V;P'?cj,cJc.C., 

pqrs pqrs 


(5.27) 


According to (5.21) the only non-vanishing contraction comes from operators of the form pp\ 
first creating a hole or a particle in the state p and subsequently destroying it. Thus, the first 
term in (5.27) must vanish. The second and the third term can only survive if p,q,r and s 
are hole indices i, j, /c and I, giving 

(4.|fli(0)|<l.) = 1 = -\T.Y + \T. Y- ('i'28) 

ij kl ij kl ij ij 


Note that neither of the two terms individually respects the Pauli exclusion principle. How¬ 
ever, the offending terms, where i = j, cancel in the sum of all terms. By the merit of Wick’s 
theorem it is no longer necessary to keep track of disallowed states individually. They simply 
cancel in the sum of all contractions. 

Let us now proceed evaluating the numerator of (5.19) in the simplified case of above 
with the perturbation given by (5.25). The next order in the expansion of the time evolution 
operator is 

0 

U^'^\0,-oo) =-i J dtiHiiti). 

— OO 


For the application of Wick’s theorem at t 7 ^ 0 it is more convenient to write the Hamiltonian 
of the perturbation Hi[t) in terms of time dependent creation and annihilation operators. 
All states in (5.25) appear in both exponents of PHTot except for p, q and r, s which can only 
appear in the left or in the right exponent, respectively. Thus, we get 
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which can be used to evaluate the next order: — oo)|<1>) 


= d> 


i f Yj ^FrKv4c\crCscl{tl)cl{tl)Cv{tl)Cnj{tl) 


pqrstuvw 


d> 


There are 4 creation and 4 annihilation operators in this expression, thus there are 4! = 24 
non-vanishing ways of contractions possible. We will, for now, just look at the following: 


—1 


f 

' —OO 

rO 




pqrstuvw 


i f Y VsF'‘Kv4>^w{ti)cscl{ti)c\cv{ti)crcl{ti) 


' —OO 


pqrstuvw 


For non-vanishing terms, p,w and q,v must be hole indices i and j, respectively. Similarly, 
s, t and r, u must be particle indices a and b, giving 


pO 1 

^ ijab 


sr 


ab tj 


4 Si -)- Sj — Sa — Sh -j- ir/ 


(5.29) 


For a system where there exists a finite gap Eg > 0 such that \£i — ea\ > Eg for all holes i and 
particles a we can simply choose rj Eg to get a result independent of the choice of rj. For 
metals with a non-degenerate ground state there is no finite > 0 for all i and a, however, 
\£i — £a\ >0. In this case the limit ?? —)■ 0 can only be taken after summing over all states, as 
done in Section 6.4 for the Uniform Electron Gas. If the ground state is degenerate we have 
to resort to Degenerate State Perturbation Theory, which is not discussed here. 

Wick’s theorem provides a systematic way to evaluate the vacuum expectation values 
occurring in the Gell-Mann-Low theorem (5.19). Instead of keeping track of which states are 
occupied after each interaction we can simply sum over contractions of the operator matrices 
Vir and Vq occurring in the perturbation, as shown in (5.28) and in (5.29). We do, however, 
have to sum over all possible contractions. 


5.4 Goldstone diagrams 

The number of possible contractions in Wick’s theorem quickly becomes too large to evaluate 
vacuum expectation values as we did in the previous section. So far, we have used diagrams 
solely for depicting the action of a second quantized operator where the initial state was 
shown below and the final state above the operator symbol, as shown in Figure 5.1. It is 
time to rigorously introduce the diagrammatic notation employed by Goldstone in order to 
enumerate and evaluate all contractions. 

Each Goulomb interaction V^e is represented by a horizontal wiggly line between two 
vertices. Each vertex consists of one electron creation operator and one electron annihilation 
operator represented by an outbound and an inbound leg, respectively. Without loss of 
generality, the outer operators cj, and Cg are associated with the left vertex and the inner 
operators Cg and Cr with the right vertex. Each effective interaction I4ff is represented by a 
vertex in form of a shaded circle. Therefore, lower indices of the matrices Vir and Vq are 
from inbound legs, upper indices are from outbound legs and left indices are from left legs. 
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Vi^CpCqCrCs — 



PA'- 

V^c'pCq = 



Figure 5.2; Representation of the Coulomb and effective interaction in Goldstone diagrams 


Connections of these legs represent the contractions, where only outbound legs may be 
connected to inbound legs, as they are the only non-vanishing operator contractions. The 
connections are directed from an outbound to an inbound leg, indicated by an arrow. Time 
propagates from bottom to top and the last Coulomb or effective interaction is usually at 

.-It 

t = 0. Connections directed forwards in time represent contractions of the form Cp.. .Cq 
and are hence restricted to particles. Conversely, connections directed backwards in time are 

restricted to holes. Connections between legs of the same Coulomb interaction are called 

R-f-' 

non-propagating and are also restricted to holes, since they are of the form CpCgCrCg. Each 
connection is labeled with a unique index, over which to contract, a,b,c,... for particles and 
i, j, k,... for holes. 

1 ^ 1 I n I 

2 2 ^ ^ OWD j 

ij j ij 

Figure 5.3: Goldstone diagrams for the terms of (5.28) with non-propagating connections 


5.4.1 Symmetries 

Interchanging the left and the right side of a Coulomb interaction including all its connections 
in the Goldstone diagram leaves the the respective matrix elements Vsr invariant. In general, 
however, we get a different contraction from the swapped connections. Since we have to sum 
over all possible contractions, both contributions have to be counted if they are distinct. Fig¬ 
ure 5.4 shows all possible left/right interchanges for the diagram representing the contraction 
evaluted in (5.29). For the sake of simplicity the time arguments of the operators are omit¬ 
ted. The two cases on the left hand side have distinct contractions, so both must be counted. 
Both evaluate to the same number since they can be transformed into each other with an 
even number of operator transpositions and Vsr = Vrs- We can do that with every Coulomb 
interaction and from now on identify one Goldstone diagram with all contractions that arise 
from interchanging the left and the right side of all its Coulomb interactions. In general, 
this cancels the factor ^ in all Coulomb operators for evaluating one Goldstone diagram as 
opposed to evaluating one single contraction, where this factor is needed according to (5.7). 

If there is however a global left/right symmetry in the Goldstone diagram, interchanging 
all Coulomb interactions simultaneously does not give a distinct contraction. It merely inter¬ 
changes the names of the indices over which to contract, as shown in the right two cases of 
Figure 5.4. Thus, for a Goldstone diagram with a global left/right symmetry only half of the 
cases arising from interchanging each Coulomb interaction are distinct, which gives rise to a 
factor of ^ for the whole diagram. 
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Figure 5.4: All possible left/right interchanges of Coulomb interactions for the diagram repre¬ 
senting the contraction evaluated in (5.29). The sums and the time arguments of the operators 
are omitted. 

5.4.2 Time integration 

The number of interactions occurring in the diagram is called its order, which is n -|- 1, where 
n is the order of the expansion of the time evolution operator in (5.17). According to 
the expansion of Ur^, we must integrate over all times ti,... ,tn of all interactions except the 
last one, respecting the order 0 > > ... > We can make the substitutions toi = 0 — 

ti,..., tn-in = tn-i—tn to integrate over the times between each interaction toiTi 2 , • • •, tn-in- 
For a particle state a propagating from an interaction at the time tn to the time t = 0 this 
gives for instance 


... e*?*! ... ... = 

0>tl>...>tn 

oo oo oo 

0 0 0 

...a...a^.. (5.30) 

So the eigenenergy Sa appears in the exponential of each time interval where the state propa¬ 
gates. For hole states i the sign is inverted. Integrating out all time intervals gives an energy 
denominator for each interval between two interactions of the form 

1 

£a + kit] 

where and Pk are the sets of the holes and particles propagating in the respective A:-th 
time interval, k is counted from the bottom. Figure 5.5 shows an example with two holes 
and two particles. 

For finite orders n of the diagrams and systems with a finite energy gap Eg one can choose 
nr] <C Eg to retrieve results independent of rj. When going to infinite orders in the diagrams, 
as it is done in the Random Phase Approximation, one has to care of the limit of ?7 —)■ 0. 
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1 

2 ^ ei+ Ej - Ea- eb + IT] 

1 onh 



{a, 6} 


Figure 5.5; Evaluation of all contractions represented by the given Goldstone diagram. There 
is one time interval where the holes i, j and the particles a, h are propagating. The symmetry 
factor is 


5.4.3 Fermion sign 

Finally, we need a rule to determine the sign of all contractions arising from one Goldstone 
diagram in a graphical way. The sign of a contraction depends on whether an even or an odd 
number P of transpositions is required to reorder the contractions in pairs: 


I I—I-1 I—I I—I 

ABC...Z = {-lfACBZ... 

Note that contracted operators must not be swapped since that turns a particle connection 
into a hole connection or vice versa, representing a different diagram. In the upper example, 

I—^ [Wl • 

the following reordering is therefore not allowed: ACZB ... In Goldstone diagrams the con¬ 
tracted operators are the operators occurring in the perturbation. They are of the general 
form 

. . . C^pC^qCrC-s ■ ■ ■ V^/Cp/Cq' ■ ■ ■ 

omiting the sums for brevity. Their respective time arguments can be integrated out according 
to the last section, the order of the interactions from right to left is, however, still relevant. 

To determine the sign of all contractions of one Goldstone diagram, we first reorder the 
operators such that operators of the same vertex of an interaction are grouped together: 

. . . VJ!^ cjjCs c]jCr ■ ■ ■ V^iCpiCq' ■ ■ ■ 
vertex 


This only affects the Goulomb interactions and involves an even number of transpositions 
leaving the sign invariant. The matrix elements Vsr and Vq are complex numbers. Their 
order is therefore irrelevant and we will omit them here. 

The finite set of vertices and connections form a directed graph, where each vertex has 
exactly one outgoing and one incoming connection. Thus, the graph consists of disconnected 
loops of vertices allowing us to group the vertices of each loop together as long as we do not 
change the order of the vertices within each loop. The sign will not be affected since we only 
move vertices. For the diagram in Figure 5.5 this gives for instance 


I n I I n J 
iajba^i^Pj^ = (+1) iaa^jbP . 


(5.31) 


The symmetry of the Coulomb interaction allows us to freely move the vertices horizontally 
such that each loop forms a simple polygon where the connections do not intersect each other. 
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Without loss of generality, we choose a clockwise orientation of the directed connections in 
order to resemble the order of the creation operators in a vertex of the form as shown in 
Figure 5.6. 



Figure 5.6; A closed loop of vertices Ai,..., Aq, where each vertex operator A is of one of 
the four possible forms shown. 


Simple polygons can be decomposed into truncated monotone polygons (Preparata and 
Shamos 2008). In a monotone polygon a horizontal line intersects the edges at most twice, such 
that its left side consists only of particle connections and its right side only of hole connections. 
In the loop in Figure 5.6 there are two monotone polygons starting at the vertices Ai and 
A 3 . After merging a particle and a hole connection from two different monotone polygons at 
vertex A 2 the rest of the polygon above the dotted line is also monotone ending at vertex A 5 . 

We have already reordered the operators to group vertices and loops together, as shown 
in (5.31). As a last preparatory step we will group vertices of monotone polygons together, 
keeping the time order within each monotone polygon. For the loop in Figure 5.6 this rear¬ 
rangement from the initial time order in the gives 

A5A4A2A6A3A1 = (+1) A5A4^A2{A3){AqAi), 

where the operators in the left and in the right parenthesis respectively belong to the left and 
right monotone polygon below the dotted line. 

With the operators ordered this way it is now sufficient to treat monotone polygons only. 
We will follow the vertices within a monotone polygon starting with a vertex of the form 
d)A. Then, vertices of the form Da add particle connections on the left side of the monotone 
polygon under consideration while vertices of the form jA add hole connections on the right 
side. A vertex of the form ia can either join two different monotone polygons, as A 2 does in 
Figure 5.6, or it can close the loop. At each vertex we reorder the involved contractions into 
the general form 


ia....a 




(5.32) 


n n 


such that P consists only of contracted pairs Wjj^... and the creation operators of the 
open particle and hole connections are to the left of P. All operators to the right of are 
then in the desired order, denoted here by a circle in the diagram. Note that the order of the 
creation operators cAA is relevant and the particle creation operator must stand to the left of 
the hole creation operator. 

A monotone polygon starts with a vertex of the form cAA. This is trivially in the above 
form 



I 


I 


ia....cAA = {+l)ia... cAA I 



(5.33) 




20 


CHAPTER 5. MANY BODY PERTURBATION THEORY 


with no change of sign. A vertex of the form b^a adds a particle connection on the left side 
of the polygon. Using two transpositions we can bring the operators into the desired form 



'' 1 n I .. I' } I n .. 

lb... b^aaU^P = {+l)ib... 



(5.34) 


such that aa^ can be absorbed into the set of paired operators P' to continue with. The sign 
does not change in this case. A vertex of the form ij^ adds a hole connection on the right side 
of the polygon. Now, three transpositions are required to bring the operators into the form 



1, 


r 


ja. 


..... . 1 i .^ . 

...ij^aU^P = (—l)ja... ii^P 



(5.35) 


to absorb ii^ into P' and to bring the hole creation operator of the open connection to 
the right side of the particle creation operator a^, as required by (5.32). In this case the sign 
changes. 

At a vertex of the form ja two different monotone polygons can be joined to a polygon 
that is monotone after this vertex. We assume that the right monotone polygon starts earlier 
in time, such that all its operators a^^P are to the right of all operators of the left monotone 
polygon Nj^Q. The converse case can be treated analogously. Bringing the operators into 
the desired form 



_■ o'' 

I hjn I I 1 1 I I n n 

ib... jab^j^QaU^P = {—l)ib... b^i^ jj^Qaa^ P 

' -V-" 

p' 



(5.36) 


requires an odd number of transpositions and changes the sign. The last vertex of a loop is 
of the form ia. It closes the monotone polygon without change of sign 


a 



.TO, . n, n ^ 

iaaH'P = (+1) aaHi'P^ . 

P' 


(5.37) 


n 

In summary, the sign changes when a contraction of hole operators ii^ is absorbed into 
the set of paired operators P', as in (5.35) and (5.36). When closing the loop in (5.37) there 
is however no change of sign although there is a pair of hole operators absorbed in P'. The 
Fermion sign of a single loop is therefore (—1)^“^ = (—1)(—1)^, where h is the number of hole 
connections in the loop. For a Goldstone diagram consisting of I loops and h hole connections 
in total the Fermion sign is thus 

(-1)'+^, (5.38) 

which can easily be determined graphically. Note that non-propagating connections, as shown 
in Figure 5.3 also count as holes. 
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5.5 The Linked-Cluster theorem 


We can now employ the framework of Goldstone diagrams to systematically evaluate the 
terms of the Gell-Mann-Low theorem (5.19): 


Eq T ^E — 


(a>|/7ot/^(0,-oo)|ch) ^ (cI>|Fi(0)17^(0,-oo)|ch) 


(<h|t/^( 0 ,-oo)|<h) (ch It/ 40 ,- 00 ) I <!>) 

Only fully contracted operators contribute to the vacuum expectation values which requires 
that all occurring diagrams are closed with no dangling connection left. In Figure 5.3 and 5.5 
we have already evaluated selected diagrams of the numerator containing the final interaction 
Hi{T). Diagrams that are connected to the final interaction at t = 0 are called connected or, 
historically, linked diagrams. In second order there are already disconnected diagrams, such 


as 


Hi{0) 0—0 


n 




dti Vki'-kWk'^ = 

ijkl 


^3 



Hi{ti) 0—0 

where the lower diagram is disconnected. Since there are no contractions between discon¬ 
nected diagrams their vacuum expectation value decomposes into factors which can be eval¬ 
uated independently. Note that each disconnected diagram comes with an additional factor 
of 1/ir/ which diverges in the limit of r/ —)■ 0. In general, we get the following form for the 
numerator of the Gell-Mann-Low energy expression containing the final interaction: 

(<I>|iLi( 0 )[/ 40 ,-oo)|d>) = (e connected diagrams^ disconnected diagrams'^ . 


) 


Let us now look at the denominator ($|f4(0, —oo)|$), where there is no operator at t = 0. 
Since the diagrams must be closed we simply get the same diagrams as the disconnected 
diagrams in the previous case, as for instance in 

t = 0 . 


Hi{ti) O- O 

Thus, the denominator has the following general form 



($ 11740 ,- 00 ) 1 ^) = E disconnected diagrams 


Finally, we need to evaluate the numerator containing Hq, which is by construction diag¬ 
onal and from (5.5) given by Hq = Eq — eiiH + Yla Therefore, all diagrams from [/^ 

must already be closed at f = 0 for non-vanishing contributions, just like in the denominator. 
In general we get 

($|i7of4(0) —oo)|$) = Eq^ Y^ disconnected diagrams^ . 

Thus, all disconnected diagrams, containing diverging factors 1 /ir/, cancel and we finally arrive 
at the expression for the correlation energy in many-body perturbation theory 


AE = 


Y2 connected diagrams j. 
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5.6 Hartree-Fock reference 

For the sake of simplicity we have not yet taken the effective interaction Ves into consideration 
although it is an integral part of the perturbation Hi. Including the effective interaction in 
first order gives just one additional diagram leading to three non-vanishing contributions: 

+ ~ (5.39) 


The perturbation Hi contains I4fr with a negative sign which has to be taken into account 
additionally to the sign of the Goldstone diagram originating from the number its loops and 
holes. The former is explicitly given here, the latter not. 

For a Hartree-Fock Hamiltonian Hq, Ves can be explicitly given 


.i.p A A 

VqCpCq 






or in diagrams: 



H-'O + 



We can insert (5.40) into (5.39) giving 


(5.40) 


In the case of the Hartree-Fock reference it turns out that the effective interaction contained 
in the perturbation exactly cancels with all Coulomb interactions containing non-propagating 
connections. This greatly simplifies the set of Goldstone diagrams to consider. In second order 
there are already 11 connected Goldstone diagrams, which are shown in Figure 5.7, and only 
two of them do not cancel in Hartree-Fock. Finite order many-body perturbation theory based 
on Hartree-Fock is referred to as Mpller-Plesset perturbation theory: MP2, MP3 or MP4. 
The maximum order is given as suffix and rarely exceeds four. Mpller-Plesset perturbation 
theory was developed by (Mpller and Plesset 1934) well before Goldstone introduced the 
diagrammatic treatment of perturbation theory discussed here. 

For a non-Hartree-Fock reference, such as Density Functional Theory, one must take all 
diagrams into consideration that contain either the effective interaction or non-propagating 
connections. There are nine such diagrams in second order, shown on the left in Figure 5.7. 
In case of Density Functional Theory, where the effective interaction consists of the Hartree 
contribution and an exchange-correlation interaction 



P 

Q 




+ 



Kc 


all diagrams containing the Hartree contribution cancel. This leaves only 4 of the 9 additional 
diagrams and they are to be evaluated with I4c instead of I4fr. Note that these diagrams 
are convergent in second order even for metals, offering an alternative to renormalization 
as done for instance by (Ren et al. 2013). Unless explicitly stated, the diagrams containing 
the effective or the exchange-correlation interaction are only taken into account to first order 
according to (5.39) computing the Hartree-Fock interaction energy with the DFT orbitals. 
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+ 










Figure 5.7: All connected Goldstone diagrams of second order. In a Hartree-Fock reference 
the left three columns cancel leaving only the two contributions of second order M0ller-Plesset 
perturbation theory (MP2). 


5.7 Propagators 

In the current approach to many-body perturbation theory we use the matrix elements Vsr 
and Vq from the second quantized representation of the interactions present in the perturbation 
Hi. However, storing the Coulomb integrals Vsr on the computer requires a large amount of 
memory scaling like O(iV^) with the size of the system since this matrix has four indices. It 
is also time consuming to calculate all elements of Vir scaling like 0{N^) with system size. 
Vir must be computed before we can even begin to evaluate any of its contractions for the 
diagrams in the perturbation expansion. 

It can be beneficial to defer the calculation of the Coulomb integrals to a later stage, espe¬ 
cially for diagrams where the sums of the Coulomb integrals can be factored into independent 
contributions, such as in the Hartree term 

O-WD = 

2 S// dxidx 2 (X 2 ) 

ij ij 


which can even be computed in 0{N‘^). We can rewrite the above expression to 


1 

2 


dxidx2 Go(xiO, xiO) ^-r Go(x20, X 2 O) , 

|ri - r2| 


defining the complex valued propagator or Green’s function from the spin-position x' to the 
spin-position x at the same instance in time t 


Go(xt, x't) ■= • 

i 


(5.41) 
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Note that the sign in the non-propagating case for equal times is negative in accordance with 
the Fermion sign rules discussed in Subsection 5.4.3. 

For propagating states we can use the time dependent creation and annihilation operators 
to express contractions over Vsr and Vq in terms of complex valued functions of two spin, 
space and time coordinates. In case of the direct MP2 diagram this gives 




L *Jill 

V’a(xi)V’a(x3)e'^“‘ 

i a 

^ V'](x2)V’i(x4)e"'''^* ^ V'fe(x2)V’6(x4)e‘^*'* . 

j b 


According to the Linked-Cluster theorem we can choose r] arbitrarily small since the terms 
diverging with 77 — 7-0 cancel. Therefore, the adiabatic switching function has no physical 
effect and rather serves to make the integrals convergent in the considered interval. Thus, we 
can as well absorb the switching function into the time dependent creation and annihilation 
operators. We get 

O ^ ■!--JJJJ 

7A:(xi)7A,(x3)e(--“+^)* V’a(xi)C(x3)e(-“+^)* 

i a 

Y, V^](x2)7A,-(x4)e(--^'+^)* J] 7/>fe(x2)7/;,*(x4)e(-'>+’?)‘, (5.42) 

j b 

'• -V----^ 

=—Go(x4t, X2O) =+Go(x 20 , X4t) 


extending the definition of the propagator^ to all cases: 


Go(xt, xt') 


' - fort<t' 

i 

+ ^ ^ otherwise. 

^ a 


(5.43) 


In other words, particle states a propagate forwards in time form t' to t, while hole states 
i propagate backwards in time from t' to t. Note that we have to use hole states i in the 
non-propagating case, where t' = t, to be consistent with dehnition (5.41). The adiabatic 
switching function is now part the propagator and its sole purpose is to make positive time 
intervals convergent for particles and negative time intervals convergent for holes. 

^Note that it is often iGo which is defined by the right hand side of (5.43). This factor is normally introduced 
such that the time evolution of the propagator is compatible to that of the Hamiltonian. We omit this factor 
here for brevity and in accordance with (Lancaster and Blundell 2014). 
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5.8 Feynman diagrams 


Eventually, it is desirable to have a diagrammatic framework that treats space and time on 
equal footing. Here, we do not need it for a relativistic treatment of the many-body system, 
however we want to be able to work in the frequency domain. This requires independent inte¬ 
grals over the whole time domain of the form dti... dtn rather than the dependent 


integrals 




dti... dtn we have from the expansion of the time evolution operator 


?7,,(0, —oo) according to (5.17): 


OO « 

n=0 \ 


dti. ..dtnHi{ti). ..Hi{tn). 

Each interaction that comes from the time evolution operator introduces a time variable to 
integrate over and a factor of (—i). In connected diagrams that applies to all interactions 
except the last one. We can, however, introduce the same for the last interaction and rewrite 
(5.42) to 

- J dtidfa fftf dxi dx2 dx3 dx4 (5(tii — ^ 

t\>H 


(-i) 




-1 


|ri - r 2 | Ts - r4| 


Go(xitl, Xsts) Go(x3t3, xHl) Go{x 2 tl,xR^) Go(x 4 t 3 , X2tl) 


Since the propagators Go only depend on time differences we no longer explicitly require the 
times to be negative and we can integrate over the whole domain instead. However, we still 
require the times to be ordered and we also need to anchor one of the interaction times for a 
convergent result. Without loss of generality, we choose to fix the last interaction at ti = 0 
by including 6 {ti). 

In the case of the MP2 direct diagram, we can now simply drop the constraints on the 
time variables taking double counting into consideration: 


-1 



dtidts 



dxi dx2 dx3 dx4 5{ti 


(-i) (-i) 


ri - r2| |r3 - r4| 

Go(xiti, X 3 t 3 ) Go(x 3 t 3 , xHi) Go(x 2 tl, X4t3) Gq(X 4t3, X 2 tl) 


In general, most permutations of the order of the time variables will lead to distinct Goldstone 
diagrams that we all want to include. This will be discussed in Subsection 5.8.2. 

Einally, we introduce time variables on every vertex to arrive at an expression where time 
and space coordinates are treated on equal footing 


(-i) 



dtidt2dt3dt4 J J J J dx2 dx3 dx4 

Go(xitl, X3f3) Go(x3t3, Xltl) Go(x 4 t 4 , X2t2) Go{x2t2, xRa) 
5{h) 


-1 


n - ^2 


<J(H - t2) I ^ - G) 


Ts - r4| 



= V'{xiti,X2t2) 


dld2 d3 d4(5(ti 


Go(l, 3) Go(3,1) H(3,4) Go(4,2) Go(2,4) H(2,1) 


(5.44) 
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defining the propagator for the Coulomb interaction 


V(xt, x't') := . ^ S(t- t') 


r — r 


(5.45) 


and the short forms Go(l,2) = Go(xiti, X 2 t 2 ) and Jdl = f^dti f dri^^^. Figure 5. 
shows the coordinate labels used in (5.44) and it is called Feynman diagram. 

1 ni,2) . 


Go(l, 



Figure 5.8: In Feynman diagrams space and time coordinates are integrated over the entire 
domain and there is no required time order of the Coulomb interactions. Each vertex is given 
a label denoting its space and time coordinates. The propagator Go represents particles or 
holes propagating between the respective vertices and the propagator V represents Coulomb 
interactions between vertices. The depicted diagram is the one evaluated in (5.44). 

5.8.1 Effective interaction 

If the effective interaction is a multiplicative effective potential (x) = ^^efr(x)V'q(x) 

we can simply use the potential Uefr(x) to evaluate Feynman diagrams containing effective 
interactions, respecting its negative sign. For a process involving a single effective interaction 
to be inserted between the two space-time coordinates 1 and 2 we get for instance 

2 

^3 = - J d3Go(2,3)ueff(3)Go(3,l). 

1 


5.8.2 Symmetries 

When dropping the constraints on the order of the interaction times we get n! permutations 
of the initial order. We now have to consider, how many of them lead to distinct Goldstone 
diagrams and in turn to distinct contractions, which we all have to count. Due to symmetries 
some permutations may lead to the same Goldstone diagram and must not be counted more 
than once. Figure 5.9 shows the case for a third order diagram, where three of the six possible 
permutations lead to distinct Goldstone diagrams. This is due to the reflection symmetry 
when swapping ti and t 2 , indicated by the dotted line. 

The product of the order of all symmetry operations on a certain Feynman diagram is 
called its symmetry factor. In the example in Figure 5.9 there is only one reflection symmetry 
with the order 2. Thus, the symmetry factor of the considered diagram is 2, which means that 
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Figure 5.9: One Feynman diagram on the left represents all Goldstone diagrams originating 
from permutations of its interaction times ti,t 2 and t^. Due to the reflection symmetry when 
swapping ti and t 2 , indicated by the dotted line, only three of the six permutations give rise 
to distinct Goldstone diagrams. The symmetry factor is 2 and the diagram must be divided 
by that upon evaluation. 

only half of all the six permutations are distinct. The symmetry factor can be determined 
graphically or computer aided, considering all permutations of the vertices. For the diagram 
in Figure 5.8 there are 4 vertices, 2 Bosonic edges B and 4 Fermionic edges F: 

B = {{1,2},{3,4}} 

F = {(1,3), (3,1), (2,4), (4,2)}. 

The Bosonic edges are undirected since swapping the left and the right side of a Coulomb 
interaction leads to the same Goldstone diagram. The Fermionic edges are directed. The 
/ 1 2 3 4 \ 

permutation 2 2 ij^® instance one of four permutations leaving the sets B 

and F unaltered; 

r(B) = {{4,3},{2,l}} = B 

t(F) = {(4,2),(2,4),(3,1),(1,3)} = F, 

thus being a symmetry operation. For this diagram, there are two reflection symmetries of 
order 2 resulting in a symmetry factor of 4. Therefore, the diagram has to be divided by 4 
as done in (5.44). 

5.8.3 Fermion sign 

The Fermion sign of a Goldstone diagram is (—l)*"*"^, where I is the number of Fermion loops 
and h is the number of hole connections in the entire diagram. In Feynman diagrams, the 
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negative sign for each hole connection is already contained in the definition of the propagator 
Go in (5.43) for the hole case. Thus, only the number of loops I still needs to be taken into 
consideration and the Fermion sign of a given Feynman diagram is 

(- 1 )'. 


Summary 

Many-body perturbation theory provides a recipe to approximate the ground state energy 
of the fully interacting Hamiltonian given the spin orbitals V’p(x) ^.nd their eigenenergies Ep 
of a single body reference, such as Hartree-Fock or DFT. It is derived from time dependent 
Rayleigh-Schrodinger perturbation theory which makes it extensive and thus applicable to 
molecules as well as solids. We can write the expansion of the perturbation series leading to 
the ground state in terms of connected diagrams and we can use either the rules of Goldstone 
diagrams or the ones of Feynman diagrams to evaluate the individual terms in the expansion. 

One Goldstone diagram of order n expands in general to 2"' different contractions orig¬ 
inating from swapping the two vertices at each Goulomb interaction. The time order of 
the interactions is fixed indicated by drawing the Goulomb interactions as parallel wiggly 
lines. Goldstone diagrams are evaluated by contracting the occurring matrix elements of the 
Goulomb integrals Vsr ■ One Feynman diagram of order n expands in general to n! different 
Goldstone diagrams arising from all permutations of the order of the interactions. The order 
of the interactions is not fixed and Coulomb interactions are normally not drawn as parallel 
lines. Feynman diagrams are evaluated by integrating the spin, space and time coordinates 
of all its vertices which are arguments to complex valued functions, the propagators, defining 
its connections. 

The symmetry of a diagram determines how many distinct contractions arise from evalu¬ 
ating a single Goldstone or Feynman diagram. In a Goldstone diagram there can only be a 
global left/right symmetry with a symmetry factor of 2, while in a Feynman diagram more 
complex symmetries are possible. In both approaches we can define diagrams with open legs 
as building blocks to be inserted in larger diagrams. In the Goldstone approach this has to 
be done in a time ordered fashion while in Feynman diagrams it cannot be done in a time 
ordered way. 

Although, many-body perturbation theory is a valuable framework for approximating 
the ground state energy it also has several drawbacks. The number of diagrams is still 
infinite and one can only evaluate a small subset. Using building blocks iteratively allows for 
evaluating an infinite number of diagrams of a certain class. This improves the results in many 
cases but there are still infinitely many diagrams neglected that can not be constructed by 
iterating building blocks. The Random Phase Approximation is a prominent example for this 
procedure. Another drawback of MBPT is that it requires not only the reference solutions 
of the unexcited states 'ipi{x) and e* but also of the theoretically infinite number of excited 
states ipaix) and Ea, called virtual orbitals. The convergence with respect to the number of 
virtual orbitals if often slow and it easilly exceeds twice the number of electrons. Finally, 
MBPT is not variational and one cannot give an upper bound for the ground state energy 
as it is given by the Hartree-Fock approximation. The non-variational nature of MBPT also 
considerably complicates the evaluation of analytic gradients, for instance with respect to the 
atom positions, i.e. forces. 



Chapter 6 


The Random Phase 


Approximation 


The Random Phase Approximation (RPA) is one of the most prominent methods beyond 
Hartree-Fock or DFT. Historically, it has been derived within two different frameworks rather 
independently. Within Rayleigh-Schrodinger perturbation theory, Heisenberg already noticed 
that certain processes are diverging in the uniform electron gas due to the vanishing band gap 
and the sign of the divergence is alternating with the order. In the diagrammatic notation 
introduced later by Feynman and Goldstone these processes are 



and they are referred to as ring diagrams. These divergencies pose serious problems as they 
render any finite order perturbation theory useless for metals. However, (Macke 1950), a 
student of Heisenberg, found a finite sum of all such diagrams - later to be termed RPA - if 
one carries out the summation over the perturbation orders before summing over the states 
in the perturbation expression of each term. This reconciles the use of perturbation theory 
for metals again. The summation over the perturbation orders can be done by iterating the 


diagram (L 


as a building block like in a geometric series. 


With the advent of Quantum Field Theory the procedure of redefining summation orders 
became more common and was termed resummation or renormalization. We will not argue 
in depth whether this procedure is justified. After all, the sum of a conditionally convergent 
series depends crucially on the order in which the individual terms are summed and - even 
worse - any result can be achieved just by choosing an appropriate order. However, we can 
argue that the notion of perturbation order is in a way arbitrary regarding that it solely 
originates from iteratively solving the equation of motion for the time evolution operator 
U{t,to) in (5.14). Therefore, it seems a natural choice to sum over the perturbation order 
first. 

In 1953, Bohm and Pines developed the RPA independenlty in the framework of the 
adiabatic connection (AC) and coined the term Random Phase Approximation. Both frame¬ 
works arrive at the same result in case of the RPA but despite the common use of diagrams 
in the adiabatic connection their meaning differs from Goldstone and Feynman diagrams of 
many-body perturbation theory discussed in Chapter 5. 
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From a more applied point of view, RPA poses an important improvement over Hartree- 
Fock and DFT. Just like finite order perturbation theories, such as Mpller-Plesset PT, it 
can describe van der Waals interaction but unlike its finite order counterparts it can also be 
applied to metals. Recent developments allow the RPA to be calculated in 0{N^) steps, just 
as DFT but with a considerably higher prefactor compared to DFT (Kaltak, Klimes, and 
Kresse 2014b). 

6.1 RPA in the frequency domain 

We will first derive the Random Phase Approximation using many-body perturbation theory 
in the frequency domain. Given the propagators Go(xt,xT') and R(xt,xT') of a system 
according to (5.43) and (5.45), we introduce the matrix notation^ 

Goxx'(^ - ^0 = G'o(xt,xY) and Vxx'(^ - t') = ^ - t') (6.1) 

I r — r I 

noting that the propagators only depend on the time difference. We define the trace and the 
matrix product of two propagators A and B by 

tr{ A} = y dx , (AB),,» = j A,,-B,-,» . 

Next, we define the independent particle polarizability xo as the first building block of the 
ring diagrams and we let Xq denote its matrix representation 

0 = xo(xt,xV) = -Go(xt,xY) Go(xY,xt), Xoxx'(^ - ^0 = Xo(xt,xY) (6.2) 

The negative sign is required according to the Fermion sign rule of Feynman diagrams since 
Xo is one closed Fermion loop. 

From now on, we will connect the building blocks Xq and V in series. In the time domain 
this corresponds to convolutions while it corresponds to simple products in the frequency 
domain. Thus, we transform Xq into the frequency domain with respect to the time difference, 
giving 

Xo(a;) = y d(t - t') e+‘^d-i') Xo(t - t'). 

We use for the forward Fourier transform into the frequency domain. Using this 

convention, the poles of the polarizability as a function of cn coincide with the positive elemen¬ 
tary excitation energies £a — £i of Hq. See (6.41) for more details in the case of the uniform 
electron gas. The Coulomb propagator V is independent of the frequency. 

We can now evaluate the diagrams of the Random Phase Approximation. In (5.44) we 
already derived an expression for the second order ring diagram in the time domain where 
di = dxj dtp. 

(-i)A A = ] Jill dld2d3d4J(ti) Go(l,3)Go(3,l) U(3,4) Go(4,2)Go(2,4) U(2,l) 

=-Xo(l,3) =-X0(4,2) 

^ A matrix with continuous indices is actually an operator. However, for numerical evaluation the coordinates 
will be discretized justifying a matrix notation for the practical application. 
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In the frequency domain there is only one frequency oj to integrate over, since we have only 
one loop and the frequency is conserved at every vertex. The integral over all frequencies 
corresponds to a convolution of all time differences rather than absolute times. Therefore, we 
do not need to anchor the diagram at a certain time anymore, here done by Using the 

matrix notation for the integrals over space we get 

= j/^tr{x„HVX„Hv} = iy'^tr{(X„{a,)V)2}. (6.3) 

This diagram has two reflection symmetries of order 2. Thus, the symmetry factor of this 
diagram is 4 giving rise to the factor 1/4 as discussed in Section 5.8. The next diagram in 
the RPA is the third order ring diagram. It has one reflection symmetry of order 2 and one 
rotational symmetry of order 3 indicated by the dotted lines. It is therefore given by 



All ring diagrams have a reflection symmetry due the symmetry of the independent particle 
polarizability Xq. Additionally, each ring diagram of order n has a rotational symmetry of 
order n, allowing us to evaluate any given order 

We can now do the resummation, summing over the orders before evaluating the frequency 
integration and the trace: 


(6.5) 



( 6 . 6 ) 


For the diagrammatic notation of the series we use the screened interaction W in RPA given 
by 

W 

H + 





+ ... 


(6.7) 


Note that the graphical appearance of a diagram containing the screened interaction W might 
be deceptive. The eye suggests a simple reflection symmetry of this diagram but in fact the 
symmetry factor must be considered for each order separately, as we have done it here. 

Instead of carrying out the matrix products in (6.6) order by order we search for a matrix 
function having the same series expansion. The function log(l — x) has the power series 
—X — x‘^/2 — x^/3 — ... so we can write the RPA energy as 



( 6 . 8 ) 
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Xo (oj) has poles along the real frequency axis making a numerical integration difficult. We 
can rotate the integration contour as long as we do not cross any poles and use the imaginary 
frequency instead. This rotation is called Wick rotation and it is discussed in more detail 
for the uniform electron gas in Section 6.4. Given Xo(ij^) in imaginary frequency we can 
substitute u = \v in (6.8) and finally get 

= \J ^tr| log (^1 - Xo(m)V^ + Xo(m)v| . (6.9) 

For the uniform electron gas (UEG) we can evaluate Xo(ij^) and V analytically and use 
(6.9) to evaluate the RPA energy for the UEG numerically. This is done in Section 6.4. For 
a molecule or a solid Xo(ii^) has to be computed from the Hartree-Fock or DFT spin-orbitals 
V'p(x). Instead of calculating Xo(t—t') in real time according to (6.2), one can already perform 
the Wick rotation in the time domain and evaluate 



Xoxx'(!''■) = -Go> 


/(ir)Gox'x(-i7') 


( 6 . 10 ) 


in imaginary time it. Note that the matrices are multiplied elementwise. The particle/hole 
propagator in imaginary time is given by 


Goxx'(i'r) = { 


{ - forr<0 

i 

+ ^ V'a(x)V’a(x06~^^“~^^^ otherwise, 


( 6 . 11 ) 


where fx is the Fermi energy. Evaluating the energies of second order Mpller-Plesset Pertur¬ 
bation Theory using the imaginary time propagators is equivalent to the Laplace transformed 
MP2 approach proposed by (Almlof 1991). To evaluate the Random Phase Approximation 
Xo (ir) needs to be Fourier transformed with respect to r to arrive at the independent par¬ 
ticle polarizability in imaginary frequency Xo(ij^) employed by (6.9). The Fourier transform 
from imaginary time to imaginary frequency, as well as the imaginary frequency integration 
in (6.9) can be done numerically on a non-equidistant grid to high accuracy with a only few 
integration points (Kaltak, Klimes, and Kresse 2014b; Kaltak, Klimes, and Kresse 2014a). To 
determine the employed quadrature frequencies and weights a function is chosen that resem¬ 
bles the RPA energy function and whose exact frequency integral is known. The proposed 
function of imaginary time is the direct MP2 term 


I [ dv 2{£a - £i) 2{eb - £j) 


1 


2 J 2-k {£a - £iY + (Sb - £j)^ £i + £j - £a - Sb 

which is the lowest order of the RPA expansion. The quadrature frequencies ixk and weights 
Wk can then be fit such that the dominant terms with a = b and i = j are best reproduced 
by the numeric integral 


1 


T £a 


-E 
2 ^ 

k 


Wk 


2{£a ^i) 

{£a - Eif + 


( 6 . 12 ) 


This fit is done for all single particle excitation energies £a—Ei and the quality of the fit depends 
on the ratio of the largest and the smallest excitation energy max(ea — e*)/min(ea — ej). For 
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Figure 6.1: Goldstone diagram of an RPA ring diagram forming a closed loop. Each Coulomb 
interaction is connected to the ring in one of four possible ways of which three are distinct. 


a non-metallic system the number of frequency points is negligible compared to the number 
of possible excitations, required in the conventional approach to calculate Xq from the Adler- 
Wiser formula (Adler 1962; Wiser 1963). The above frequency grid allows the RPA energy to 
be evaluated in 0{N^) steps, just like DFT but with a considerably higher prefactor. For a 
metallic system the behavior of the RPA energy for large imaginary frequencies is important 
for an accurate numerical quadrature. This is discussed in the end of Section 6.4. 

6.2 Direct Ring Coupled Cluster Doubles 

The Random Phase Approximation can also be evaluated using the matrix elements of the 
Coulomb integral 

VFr = jj dxdx>p(x)V’*(x') 1^:3^ V'r(x')V's(x) 

rather than the propagators Gq and V. This approach is not as efficient, however an important 
correction to the error remaining in the RPA is based on this approach. 

The RPA ring diagrams form a closed loop and Figure 6.1 shows the Coldstone diagram of 
one such ring diagram. As discussed in Section 5.4, connected vertices represent contractions 
and we need to integrate over the time interval between Coulomb interactions. We will do 
that from bottom to top following the left and the right particle/hole pairs along the loop. 
In the diagram given as example in Figure 6.1 we start with interaction 1 following the left 
and the right particle/hole pairs to interaction 3. There, the right pair is contracted at one 
vertex and a new pair to follow on the right side emerges on the other vertex of interaction 3. 
Independently, we also follow the left and the right particle/hole pairs starting at interaction 
2. At interaction 4 the two processes merge such that there is only one left pair and one right 
pair remaining after the time indicated by the dotted line. At interaction 5 both pairs are 
finally contracted and the loop is closed. This procedure is analogous to the one employed 
for deriving the Fermion sign of Coldstone diagrams in Subsection 5.4.3. 

While following the left and the right particle/hole pairs we perform all occurring con¬ 
tractions and integrals over the respective time intervals and we keep the intermediate results 
in a matrix depending on the states of two particle/hole pairs, diagrammatically denoted 

aiiJ 

= v.v ■ 

is called direct ring Coupled Cluster Doubles (drCCD) amplitudes and it contains the 
probability amplitude to arrive at two particle/hole pairs in the given states in any of the 
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possible ways forming the ring diagrams. We will now go through all cases that can occur 
following the left and the right particle/hole pair from bottom to top indicated in Figure 6.1. 

In the first case, a Coulomb interaction creates two new particle/hole pairs as in interaction 
1 in Figure 6.1. This can occur at any time in the past so we need to integrate over the time 
interval between the Coulomb interaction and the time where we want to use the probability 
amplitudes tfj*, which we always move at t = 0. One contribution to the drCCD amplitudes 
is thus 

a i b J a i b J 


t = 0 


V-.y 


t» = (-i) 




ab 

^3 


+ .. 





Si + Sj — Sa — Sh + '\rj 


+ ... 


(6.13) 


In the next case, a Coulomb interaction contracts the right particle/hole pair creating a new 
one. This is the case for interaction 3 in Figure 6.1. We can now use the drCCD amplitudes 
recursively containing all processes until the time of the interaction. The time interval between 
the interaction and the time where we want to use the new drCCD amplitudes still needs to 
be integrated as before. The second contribution to the drCCD amplitudes is thus 
a i b 3 a i b 3 


t = 0 


v__y 


j.ah _ 

^ij — 


+ 


Ekct^v. 


kb 

cj 


Si T Sj Sq 


Sb + i?? 


+ ... 



The Fermion sign of this contribution is positive since there is one more hole and one more 
closed Fermion loop. Note however, that the denominator is negative giving rise to the 
alternating sign when evaluating the RPA order by order. The same case can occur for the 
left particle/hole pair giving 


t = 0 


a i b 3 

v__y 


• • • W 


E -\rakj.cb 
kc ^ic ^kj 


Si + Sj — Sa — Sb + il] 


+ 



The last case occurs when two independent processes merge at a Coulomb interaction, as in 
interaction 4 of Figure 6.1. In this contribution, the drCCD amplitudes occur in a quadratic 
form on the right hand side 

a i b 3 a i b 3 


t = 0 


V__V 


j.a6 _ 


E ^acf/klfdb 
klcd ^ik ^cd ^Ij 

Si + Sj — Sa — St IT] 



t = 0 

(6.16) 


Taking all contribution (6.13) to (6.16) into consideration and taking the limit ry —)■ 0 yields 
the drCCD amplitudes equation 


{si + Sj - Sa- Sb) tij — + 




kc 


jrak^cb 
^ic ^kj 


kc 


klcd 


acjrkl^db 
T] 


cd ^Ij 


(6.17) 
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This equation is quadratic and can only be solved by iteration. The convergence with respect 
to the number of iterations is, however, fast. Employing a Shanks transform, 8 iterations are 
sufficient to yield a converged RPA energy to 5 significant digits of precision even for a system 
with a low band gap, such as a finite size uniform electron gas (Freeman 1977; Shanks 1955). 
Each iteration is still costly requiring 0{N^) steps. The iteration process also requires the 
amplitudes to be stored demanding 0{N^) of memory. 

Given the drCCD amplitudes, we can evaluate the RPA energy by contracting both par¬ 
ticle/hole pairs with the last Coulomb interaction, corresponding to interaction 5 in Figure 
6.1. The drCCD amplitudes have a left/right reflection symmetry. The RPA energy is thus 



where the two closing Fermion loops result in a positive Fermion sign. 

In this approach the time order is always maintained by the way the drCCD amplitudes 
are recursively defined. Therefore, no particular symmetries have to be considered apart 
from the reflection symmetry. Ring diagrams that have more than two particle/hole pairs 
at some instances in time arise from the quadratic contribution (6.16). This is the case for 
the diagram in Figure 6.1 between interaction 2 and 4. Excluding the quadratic contribution 
(6.16) gives the Tamm-Dancoff Approximation (TDA), which is the subset of all RPA ring 
diagrams where there are exactly two particle/hole pairs at all times between the first and 
the last interaction. From (6.16) and (6.13) follows that the lowest order diagram of RPA 
that is not part of TDA is of fourth order. 

On the other hand, the drCCD amplitude equations are a subset of the amplitude equa¬ 
tions including all possible ways to arrive at two particle/hole pairs. These amplitude equa¬ 
tions are called Coupled Cluster Singles Doubles (CCSD) amplitudes and they contain pro¬ 
cesses such as 



Unfortunately, it is computationally more time consuming to calculate the full CCSD am¬ 
plitudes, scaling like 0{N^), since the employed tensor multiplications cannot be split into 
pieces, involving no more than 5 indices. For the calculation of the drCCD amplitudes this 
can be done, since the matrix of the Coulomb interaction Vlr can be decomposed into a 
product of two tensors with 3 indices in the momentum basis: 

yj? = j X?(G) Xq*{G) , with xP(G) = y dx V’p(x)e'^'^V'g(x) • 

This allows tensor products involving Vir to be “cut” at the Coulomb line, accelerating the 
evaluation of the direct ring Coupled Cluster Doubles amplitudes to 0{N^). For the first 
non-trivial case in (6.17), this is done, for instance, by 

'^i“(G) = ^ t“fcXc(G) , {Si + Ej - Ea- Eh) t^j = . .. + J (G) + • • • 

The Coupled Cluster method was developed by (Coester and Kiimmel 1960) for the atomic 
nucleus and later adopted for electronic correlation by (Cfzek 1969). 
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6.3 RPA from the Adiabatic Connection 

In the framework of many-body perturbation theory the derivation of the Random Phase 
Approximation is straight forward. According to the last two sections it can be done either in 
the frequency domain employing Feynman diagrams and regarding the symmetries of the ring 
diagrams or in the time domain employing Goldstone diagrams. In Chapter 5 we discussed 
the equivalence of the two approaches. 

Despite the straight forward derivations and the general applicability of many-body per¬ 
turbation theory to arbitrary reference systems, neither of the two previously discussed deriva¬ 
tions are considered standard approaches to the Random Phase Approximation. In this sec¬ 
tion we discuss the framework of the Adiabatic Connection (AC) employed by Bohm and 
Pines, which is considered the standard approach to the RPA, at least in solid state physics. 

In the Adiabatic Connection (AC) we define a Hamiltonian H (A) depending on a coupling 
constant A specifying the strength of the full electron-electron Coulomb interaction 

^(A) =f+ I4e + I4fr(A) + AI4e, Ag[0,1]. (6.19) 

The effective interaction i4ff(A) also depends on the coupling constant and we choose I4fj(A) 
such that the density of the system with the Hamiltonian H (A) is for all A equivalent to the 
density of the fully interacting system with the Hamiltonian H{1). This is in contrast to 
many-body perturbation theory where the effective interaction is simply scaled by the factor 
1 — A according to (5.15) and where A is time dependent. Requiring a constant density for all 
A is a strong condition implying that the density of the reference system with the Hamiltonian 
H{0) is equivalent to the density of the fully interacting system. Although this condition is 
met to a good degree by using a DFT Hamiltonian for H{0), the Hohenberg-Kohn theorem 
only states the existence of such an effective potential I4fr(A). For practical considerations it 
is known that in general no DFT density fully agrees with the density of the respective fully 
interacting system. 

Let now 'I'(A) denote the normalized ground state of the respective Hamiltonian H{X), 
being the solution of 

iL(A)|T(A))=E|T(A)) . 

The ground state energy E is equivalent for all A since E isa functional of the density according 
to the Hohenberg-Kohn theorem and the density is the same for all A. Even assuming that the 
density n(r) of the DFT reference is exact we can only directly evaluate the nuclei-electron 
potential energy E^e as the functionals for the kinetic energy T and for the electron-electron 
potential energy E^e are unknown: 

E[n{r)] = T[n(r)] Ene[n{r)] + Eee[n(r)] . 

However, we are only interested in the sum of the kinetic and the potential energy so we may 
as well use the known kinetic energy of the Kohn-Sham system 

^ V2 

rs = -(T(o)| J;^|t(o)) 

n=l 

and ask for the energy to be added to Ts and Ene to arrive at the same total energy E. This 
energy is called Hartree-exchange-correlation energy and it is given by 


E-Rxc = E - Eae - Ts . 
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Using (6.19) we can express the terms above by expectation values of the fully and the non¬ 
interacting Hamiltonian: 


H - Hne = (^(1)|H(1)|^(1)) - (^'(l)|U„e + Ueff(1)|^(1)) 

Ts = ('I/(0)|H(0)|'I'(0)) - {^mVne + Uefr(0)|'I^(0)), 

noting that t4fr(l) = 0 and t4fr(0) is the effective potential of the Kohn-Sham system. This 
allows us to make the connection between the fully interacting and the non-interacting system 


dA ^ - Ke - Ueff(A)|'l/(A)) 


( 6 . 20 ) 


'I'(A) is the ground state of H{X). As a consequence, we can use the Giittinger theorem 
(Giittinger 1932), also known as Hellman-Feynman theorem, to evaluate the total derivative 
of the expectation values involving the Hamiltonian: 


^ (T(A)|H(A)|T(A)) = (t(A)| A^(A)|t(A)) 


( 6 . 21 ) 


Furthermore, l^e and 14fr(^) are local potentials and their expectation value only depends 
on the electron density n(r). Since the density is constant for all A the total derivate of their 
respective expectation value simplifies to 


^ (T(A)|Ueff(A)|^(A)) = (t(A)| A i4g(A)|T(A) 


( 6 . 22 ) 

(6.23) 


Inserting the total derivatives into (6.20) yields the final expression for the Hartree-exchange- 
correlation energy in the adiabatic connection: 


^Hxc= / dA(T(A)|Uee|^'(A)). 

Jo 


(6.24) 


In other words, the Hartree-exchange-correlation energy is the coupling strength averaged 
potential electron-electron energy. 


6.3.1 Fluctuation dissipation theorem 

Equation (6.24) cannot be evaluated directly since the ground states 'I'(A) are not available 
for any A except 0. However, the whole information of the many-body wavefunction T(A) is 
not required to evaluate the Goulomb operator. The pair density n^(x, x') is sufficient: 


(T|Uee|^) = / dxi. 


dxArT*(xi,...,XAr)^ ^ ^ 


n^m 


Tr,. - Tr. 


T(xi,... ,XAr) 


= 4 y dxy dx' ^ ^ y dxi... y dxAr T*(. ..) ^ ( 5 (x - x„)(i(x' - Xm) ^'(. • •) 

n^m 


= n2(x,x') 
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The pair density can be written in terms of unrestricted summations 


n 




'(x,x') = ^ (5(x - x„)(5(x' - x„ 

n^m 

^1 “ X„)5(x' - Xm) ^ (5(x - x')(5(x - x„ 


'If 


= (5(x—x')n(x) 

which allows us to relate it to the density fluctuation operator 5h{x) = Xln 


since 


(T|(5h(x)(5h(x')|'i;') = (|'h| ^5(x - x„) ^(5(x' - x^) + n(x)n(x') 

n m 

- ^(5(x -x„)|t^ n(x') - n(x)n(x'). 


= n(x) 


The pair density and the density fluctuation operator are thus related by 

n^(x,x') = (T|(5h(x)(5h(x')|T) + n(x)n(x') — 5(x — x')n(x). (6.25) 

Note that although ('I'|(5h(x)|T) is zero, the expectation value of a quadratic form of the 
density fluctuation operator is in general non-zero. The fluctuation dissipation theorem links 
the density-density fluctuations at the two sites x and x' occurring in (6.25) to the density- 
density response function XA(x,x',iiz) of the system with the coupling strength A, here given 
in imaginary frequency: 

-/ ^XA(x,x',ii/) = (T(A)|(5h(x)(5h(x')|T(A)). (6.26) 

With time dependent density functional theory (TDDFT) we can finally connect the density- 
density response function xa of the system with a coupling strength A to the density-density 
response function xo of the DFT reference system 


XA(xi,X4,m)= xo(xi,X4,m) 

+ [[ dx2dx3Xo(xi,X2,m) 


A 


|r2 - rsl 


+ /4(x2,X3,m) XA(x3,X4,ii/), (6.27) 


where the Coulomb interaction is scaled by A. is called exchange-correlation kernel and 
it contains the change of the DFT effective potential with respect to a change in the electron 
density, is not explicitly known and it is coupling strength dependent. (6.27) is an implicit 
integral equation for xa called Dyson-like equation. It requires that the density response of 
the interacting system to a change of the external potential is the same as the density response 
of the DFT reference system to a change of the effective potential. This assumption is the 
TDDFT counterpart of the assumption that the DFT density is exact. 

The response function of the DFT reference system xo can be evaluated from the DFT 
spin-orbitals '0p(x) with first order perturbation theory: 


Xo(x,x',m) = 


-E 


V’a(x)V’a(x')V’i(x')V’r(x) , V'*(x)V’i (x')V'a(x')V'a(x) 


Ea — £i — m 


£a- ei + \v 


(6.28) 
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See (6.45) for an analogous derivation of the above result for the uniform electron gas. Note 
that the definition of xo given in this section expects the Coulomb propagator to be l/|r —r'|. 
This differs from the convention used elsewhere in this work and has been chosen for consis¬ 
tency of this section with literature. In (6.51) this is discussed in detail. The above definition 
further assumes a spin-unrestricted reference and an imaginary frequency integration over the 
entire domain from —oo to oo. 

Using (6.24), (6.25) and (6.26) the Hartree-exchange-correlation energy is given by 


^’Hxc = 



dx 


/ (^(x)n(x') - S{x - x')n(x) - j ^XA(x,x',m) 


= En + 


1 

2 



dx 




du 


X\i 


X, X ,11/, 


(6.29) 


where the Hartree energy En is readily factored out. To tackle the remaining Dirac delta 
function we rewrite n(x) = Xli V’i(x)V'i (x) in the context of the delta function and then use 
the completeness of the eigenstates V'p(x): 


(5(x-x')n(x) =(I(x-x')^V'i(x)V’i (x) = ^ V’p(x)V’p(x')V'i(x')V'i (x) 

i pi 

= X] V’i(x)V’j(^OV'i(x )V'i(x) -h J]]V'a(x)C(x')V’i(x')V'i(x) . (6.30) 

ji ai 

The first sum over unexcited states can be inserted into (6.29) and factored out giving the 
exact exchange energy 


E^= — 
2 




V’i(x)V’i(x')V'i(x')V'r(x) 


r — r 


The exact exchange energy corresponds to the given first order diagram in many-body per¬ 
turbation theory with the Kohn-Sham spin-orbitals '0p(x) as reference system. 

The second sum in (6.30) over excited and unexcited states is the integral of xo(x, x',iz/) 
given in (6.28) over all imaginary frequencies. Note that there is only one pole at i/ = 
=Fi(ea ~ £i) requiring a convergence factor of with 0 < |r| <C 1 for a contour enclosing 

the pole with vanishing contribution at infinity. We choose to take the limit r —)■ O"*" from 
above resulting in a clockwise contour. The sum over the residue gives 


ia 


The above result can also be obtained evaluating xo(x, x',ir) using imaginary time propaga¬ 
tors according to (6.10) and (6.11) and then taking the limit r —)■ O'*'. Finally, we can insert 
the exact exchange energy and the above expression into (6.29) to arrive at an expression for 
the correlation energy; 

Ec = dX j dx j J ^ (xA(x,x',m) - xo(x, x, iz/)) , (6.31) 


where F^Hxc = Eu + E^ + Ec. 
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6.3.2 Random Phase Approximation of the polarizability 

The Dyson-like equation (6.27) for the polarizability Xa(x, x', ii/) can be written in the matrix 
notation introduced in Section 6.1 

XA(m) = Xo(m) + Xo(m) (aV + XxM . (6.32) 

In the Random Phase Approximation one considers only the Hartree term in the time de¬ 
pendent density functional theory derivation of the above Dyson-like equation, neglecting the 
exchange correlation kernel The polarizability can then be expanded to 

XA(m) = Xo(m) + AXo(m)VXo(m) + A2Xo(m)VXo(m)VXo(m) + ... 

Writing (6.31) in the matrix notation and inserting the above expansion gives 

Ec= -^^'dA j ^tr{XA(m)V-Xo(m)V} 

= -^^'dA j ^tr{A(Xo(m)V)' + A2(Xo(m)V)^ + ...} . (6.33) 

The explicitly known A dependency allows us to integrate A analytically, which leads to the 
final expression of the correlation energy in the Random Phase Approximation: 

-I I ^tr|i(Xo(m)V)V^(Xo(m)V)^ + ...| 

= +11^ tr{ log (l - Xo(m)v) + Xo(m)v} . (6.34) 

This result is formally equivalent to (6.9) although two very different approaches have 
been used in the respective derivations. The prescribed Hartree En and exact-exchange term 
Ex also correspond to the two first order diagrams in many-body perturbation theory using 
a non-Hartree-Fock reference. Note however, that the factors in the series expansion 

^(Xo(m)V)' + i(Xo(m)V)' + ... 

have entirely different reasons in the two approaches. In the approach using many-body 
perturbation theory they originate from the rotational symmetry of the ring diagrams. In 
the Adiabatic Connection, on the other side, they stem from averaging the potential energy 
over all A to arrive at the total energy. Thus, the two approaches can very well disagree when 
different classes of diagrams are considered. 


6.4 RPA for the uniform electron gas 


In this section we apply the framework of many-body perturbation theory as discussed in 
Chapter 5 to the uniform electron gas (UEG) and and use resulting propagators to calculate 
the Random Phase Approximation for this system. 

We choose the free Hamiltonian as a reference only containing the kinetic energy 


N 

^o = -E 

n=l 
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for a system with N electrons in a cubic box of length a and volume Ll = a^. The solutions of 
the non-interacting Hamiltonian are plane waves commensurate with the box. The normalized 
spin orbitals and eigenenergies are 


1 -i. 27r o 

V'(Tk(ar) = 5aa—;=e , eo-k = -pr , where k G —Z , and u, a G {t)i} • (6.35) 

V H 2 a 

The uniform electron gas is the limit of this system taking its volume H to infinity while 
keeping the electron density N/Ll fixed. The average volume per electron is usually given by 
specifying the the radius of a sphere of equal volume 


Q/N 


3 


(6.36) 


Vs is called Wigner-Seitz radius. In the non-interacting ground state the lowest N spin- 
orbitals are occupied and the wave numbers k of these states lie inside a sphere called Fermi 
sphere with radius /cp: 

E =^- (6.37) 

(T,|k|<fcp 

The density of the wave vectors k G ja increases with increasing box volume so we can 

approximate the above sum by an appropriate integral 


E = E 

(T,k (y 


Ll 


(2vr) 


dk 


(6.38) 


and use it together with (6.37) and (6.36) to find fcp as a function of the Wigner-Seitz radius 
r*: 

/ 'V Qtt 1 

(6.39) 




> Dtt 1 

\ki = 


^ 2E. r3- 

In the non-spin polarized case = 2 and in the spin polarized case = 1. 


6.4.1 Propagators 

We can now evaluate the propagators for the uniform electron gas according to (5.43): 


Go(art, olr'-d) 


— E r')e^ t) for t < t' 

(T,|k|<A:p 

-1- E ^ otherwise. 

. (T,|k|>A:p 


Inserting the spin orbitals and eigenenergies from (6.35) and approximating the sum over 
states by the appropriate integral from (6.38) gives 


Go(art, a'r'f') = < 


|k'|<A:p 


+ 


dk^ 

(27r)3 

dk' 


fort<t' 


I, |k'|>fcp 


[ otherwise, 

J ( 27 r )3 
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which we can transform into momentum space k with respect to (r — r'). For the case t < t' 
this yields 


— / d(r — r )e 


/N„-ik-(r-r') 


|k'|<fcF 


X , „ik'-(r-r') (-ik'V2+,,)(i-t') 


dk' 


{27rf (5(k' - k) 


J (27r)3 

|k'|<fcF 

= - 0(fcF - |k|) ^ 

where 6{x) is the Heaviside step function. In the general case this gives 

-9{kF - |k|) e(-‘kV2+r;)(i-i') ^ 

+0(|k| - kp) e(-‘kV2-r?)(i-i') otherwise, 


Go(k,t-t') = 


(6.40) 


as a function of k and the time difference {t — t'), omitting the electron spins. Finally, we can 
transform the propagator into frequency space uj with respect to {t — t'). For the case t < t' 
this gives 

- r d{t - t') e{kF - |k|) 

J —oo 

(1-O)0(A:f - |k|) i0(A:F-|k|) 


i(a; — k2/2) + r/ uj — k2/2 — irj ' 

Note that, unlike for the momentum, we use for the forward Fourier transform into 

frequency space so that the poles of the propagator coincide with the positive eigenenergies. 
The particle/hole propagator is now quite compact, even for the general case: 


Go(k,a;) = 


UJ - k2/2 + ir/k ’ 


where % = 


—rj for |k| < ^F 
+r] otherwise. 


(6.41) 


In the space and time domain we have to integrate over all spins, space and time coordinates 
of each vertex: 


dr dt. 


In the momentum and frequency domain we integrate over all spins, momenta and frequencies 
of all propagators connecting the vertices: 

Next, we also want to transform the propagator for the Coulomb interaction 

V (art, art') = -j- -r 5(t — t') 

|r — r I 

in the momentum and frequency domain. The transform with respect to {t — t') is straight 
forward and simply gives a frequency independent propagator 


V{r — r' ,uj) = 


—1 


r — r 
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Figure 6.2: Go{u}) has one pole above the real axis for |k| < kp or below the real axis for 
|k| > kp. We can smoothly rotate the integration contour from the real axis to a contour 
parallel to the imaginary axis without crossing any pole. The arcs indicate the contour at 
infinity. 

omitting the spin coordinates. The transform into momentum space does not converge, 
however, the propagator for the Yukawa potential can be expressed by a Fourier integral for 
any parameter m > 0: 


-ig-m\r-r'\ 


r — r 


(2vr) 


dqe^'i- 


,iq-(r—r') 


—dvri 


n 


(2vr) 


dciF^- 


iq-(r—r') 


—dvri 


n(q2 + m^) 


Taking the limit m —0 gives the propagator for the Coulomb interaction in momentum 
space: 


Y(q,a;) = 


dvri 

riq^ 


(6.43) 


Note that we use k for particle/hole propagators while we use q for Coulomb propagators. The 
spin coordinates remain omitted for brevity as the behavior of spin is intuitive in the UEG. 
A particle or a hole does not change spin during propagation while the Coulomb interaction 
mediates between particles or holes irrespective of the their spin. 


6.4.2 Imaginary frequencies 

Go has one pole on the positive real frequency axis. The pole is infinitesimally below or above 
the real axis depending on whether it is a particle, having |k| > kp, or a hole. On which side 
of the real axis the pole lies determines whether there is forwards or backwards propagation 
in time. Figure 6.2 shows the pole of Go(k, cu) as a function of cu with k as parameter. The 
poles of the propagators make a numeric frequency integration along the real axis difficult. 
We can, however, smoothly deform the integration contour as long as we do not cross any 
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poles. As indicated in Figure 6.2 we can rotate the contour Ci counterclockwise around the 
Fermi energy = /cp/2 to arrive at the contour € 2 - Substituting lo = fi + iu we get 


/ duj = <p dio = (p du = idv ^ 

J — 00 J J J —00 

Cl C 2 

where we assume that the contribution of a contour at infinity vanishes. This rotation of 
the integral contour is called Wick rotation to imaginary frequency. When integrating in 
imaginary frequency the side of the poles relative to the contour no longer depend on rj and 
we can as well take the limit 0 before evaluating the imaginary frequency integrals. The 
particle/hole propagator in imaginary frequency is then given by 


Go(k,ii/) = 


1 


in — /2 — fi) z/+ iAffk 


where Ae^ = k^/2 — /cp/2 . 


(6.44) 


6.4.3 Polarizability 


Next, we can evaluate the independent particle polarizability. The polarizability diagram 
has two open vertices. Let q and n denote the momentum and the imaginary frequency 
incident at the lower vertex shown in the diagram below. From the momentum and frequency 
conservation at every vertex follows that there is one pair of momentum k and frequency rj to 
integrate over and that the outgoing momentum and frequency is equal to the incoming ones. 
In the context of this diagram part k and rj are called internal momentum and frequency, 
respectively, while q and n are connected to other diagram parts and are therefore referred 
to as external momentum and frequency, respectively. We evaluate the diagram part by 
integrating over all internal momenta and frequencies starting with the analytic ij integration: 


c[,n 





Go(k + q, iry + iu) Go(k, ir/) 

ZTT 


1 


1 


+ 


27r rf + n + iAek+q g + iAck 
k + q| — kF)6{kF — |k|) 

^k+q 

6>(|k| - kF)e{kF - |k + q|) \ 
ek - ek+q + W ) 


(6.45) 


The conditions on k expressed by the Heaviside theta functions arise from the necessity to have 
one pole on either side of the integration contour for a non-vanishing result. The conditions 
in (6.45) can be unihed by substituting —k' = k + q in the second sum. This gives 


To(q,iz^) = -i 


r Hdk 

^ J ( 27 r )3 

|k|<fcF<|k+q| 


1 


£k+q £^k 


+ 


1 


- - in 


£k+q — £k + in 


(6.46) 


With effort this expression can also be integrated analytically with respect to the momentum 
k and expressed in terms of a real valued function R{q, u) where the Fermi momentum and 
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the Fermi energy are normalized to 1 and 1/2, respectively (Ziesche 2010): 


xo(q,ii^) = 


7r2 V/cp’lql^F 


2R{q, u) =1 — u ^arctan ^ ^ Cretan 


l-g/2 


u 


1 + u^ - g^/4 , 

+- . log 

2q 


+ {q/2 + 1)2 
^2 + {q/2 - 1)2 


(6.47) 

(6.48) 

(6.49) 


Note that the imaginary units are often traded between the Coulomb kernel 17(q) and the 
polarizability Xo{h ’^^) since only their product occurs in the expression for the RPA corre¬ 
lation energy. Throughout this work, except Section 6.3, imaginary units are used exactly as 
they emerge from Wick rotated imaginary frequency integrations: 


P(q) = 

dvri 
n q2 

x(q, ii^) 

II 

1 

used in this work. 

(6.50) 

P(q) = 

dvr 

Ll q2 

x(q, ii^) 

= -E(-) 

also common in literature. 

(6.51) 


6.4.4 Correlation energy 

Finally, we can evaluate the correlation energy in the Random Phase Approximation for the 
uniform electron gas according to (6.9). In the space domain we needed to convolve the 
propagators connected in series which we denoted by the matrix product. In the momentum 
domain of a homogeneous system this simplifies to a product. The RPA energy per electron 
in the uniform electron gas is thus 

=/j\j'T {log(i -xo(q,ii')r(q)) + xo(qi>')r(q)} 

= ^{ioi!(i-xo(<i,iOr(<,)) + xo(<i,iOr(<i)} (6.52) 

Note that the sum over all momenta q of the Coulomb propagator does not include a sum over 
spins. The given integral can been evaluated using a Gauss-Kronrod rule first in q then in 
v to yield the correlation energy in the Random Phase Approximation to 5 significant digits 
of precision for the non-spin polarized, paramagnetic case with = 2 and for the spin 
polarized, ferromagnetic case with = 1. The results are listed in Table 6.1 for different 
densities and they include Quantum Monte Carlo results for comparison. Although the RPA 
systematically overestimates the correlation energy it lacks only about 1 /3 of the correlation 
energy at r* = 1. (Gell-Mann and Brueckner 1957) have shown that in the limit of —)■ 0 the 
entire correlation energy is contained in the RPA ring diagrams. For low densities, however, 
the relative error becomes larger and other diagrams become important. 


6.4.5 Large momentum behavior 

In case of the uniform electron gas the polarizability can be evaluated for an arbitrary mag¬ 
nitude of the momentum transfer q allowing for an accurate numerical integration. For solids 
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paramagnetic 

ferromagnetic 

rs 

[a.u.] 

[vciEh iV] ± 

Ec 

[uiEhN] 

Ts 

[a.u.] 

[vnEh A^j ± 

Ec 

[uiEhN] 

1 

-78.799 

0.001 

-59.632 

1 

-51.893 

0.002 

-31.701 

2 

-61.801 

0.001 

-45.091 

2 

-42.416 

0.001 

-24.090 

3 

-52.759 

<0.001 

-37.214 

3 

-37.179 

0.001 

-20.048 

4 

-46.806 

<0.001 

-32.054 

4 

-33.633 

<0.001 

-17.415 

5 

-42.470 

<0.001 

-28.339 

5 

-30.992 

<0.001 

-15.520 

6 

-39.117 

<0.001 

-25.504 

6 

-28.911 

<0.001 

-14.071 

7 

-36.418 

<0.001 

-23.253 

7 

-27.209 

<0.001 

-12.916 

8 

-34.182 

<0.001 

-21.414 

8 

-25.778 

<0.001 

-11.969 

9 

-32.289 

<0.001 

-19.876 

9 

-24.551 

<0.001 

-11.174 

10 

-30.658 

<0.001 

-18.568 

10 

-23.482 

<0.001 

-10.495 

12 

-27.975 

<0.001 

-16.454 

12 

-21.698 

<0.001 

-9.391 

15 

-24.929 

<0.001 

-14.119 

15 

-19.629 

<0.001 

-8.160 

20 

-21.381 

<0.001 

-11.497 

20 

-17.156 

<0.001 

-6.758 

30 

-17.068 

<0.001 

-8.486 

30 

-14.044 

<0.001 

-5.112 

40 

-14.463 

<0.001 

-6.778 

40 

-12.099 

<0.001 

-4.156 

50 

-12.680 

<0.001 

-5.666 

50 

-10.736 

<0.001 

-3.521 


Table 6.1: Correlation energy in the Random Phase Approximation for the uniform elec¬ 
tron gas at different densities compared to Quantum Monte-Carlo (QMC) results obtained 
by (Ceperley and Alder 1980) fitted by (Perdew and Zunger 1981). For RPA energies at 
intermediate spin polarizations see (Vosko, Wilk, and Nusair 1980). 




6.4. RPA FOR THE UNIFORM ELECTRON GAS 


47 


or molecules the finite resolution used for finding the Hartree-Fock or DFT oribitals imposes 
however an upper limit Gmax up to where xo can be evaluated. Therefore, we need to know 
the asymptotic behavior of the RPA correlation energy for sufficiently large values of Gmax iu 
order to extrapolate numerical results to Gmax —^ co. Assuming that the system is sufficiently 
homogeneous at the resolution corresponding to a given Gmax, the asymptotic behavior of the 
RPA in a solid or in a molecule is the same as in the uniform electron gas. The latter can be 
derived analytically and will be outlined here. 

Given the expression for the independent particle polarizability xo of the uniform electron 
gas from (6.46) 


Xo(q,iJ^) = -i X] 


|k|<fcF<|k+q| 


■E 


f^dk 

(27r)^ 


f2dk 


1 


^k+q 

A£k,q 


+ 


1 


£k+q — £k + il' 


|k|<fcF<|k+q| 


(27r)3 Ae2 q + ’ 


(6.53) 


where A£k,q = £k+q ~ ^k = k • q + q^/2, we can trivially integrate out k for large enough 
magnitudes of q since |k| < fcp <C |q|. This yields 


xoiqui^) ~ 


-i /cf 


(?V2) 


2 I 2 

+ 


(6.54) 


Next, we can insert this approximation into the RPA energy expression (6.52) for a given, 
large q and integrate out the imaginary frequency u, getting 


du 


{ log (^1 - xo{q, (q)^ + Xo{q, iz^)k"(g)} 


(gV2) ^(g2/2)2 + l- (gV2)' - 1/2 


—-5-4-T7T + G ( —TJ 


q6 qio 


P4 ’ 


where we expanded in the variable u = 1/q at u = 0. The leading order term can finally be 
used to estimate the missing RPA energy per electron AFl/^^^(Gmax) when truncating the 
momentum integration at a finite momentum Gmax^ 


AEf^^{G 


PCX 

JGn 


dqq 



Gl 


+ 0 



(6.55) 


6.4.6 Large imaginary frequency behavior 

Although most implementations of the RPA can evaluate the independent particle polariz¬ 
ability Xo at arbitrary imaginary frequencies, knowledge of the asymptotic behavior of the 
RPA energy for large frequencies u is useful for choosing an appropriate variable transform 
for integrating the tail. This is relevant for metallic systems, where we can assume that the 
system behaves like a Uniform Electron Gas at times short enough. 

Unlike in the case of large momenta q, k cannot be trivially integrated out. However, we 
can separate the momentum integration of q into three regions and show that they all come 
to an analogous form of the integrand for that momentum integration. In Subsection 6.4.5 
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Figure 6.3; Cross section of the set of momenta k such that |k| < fcp < |k + q| for small 
magnitudes of the excitation momentum q. In this limit, the set covers half the surface of 
the Fermi sphere and its maximum thickness is |q|. 


we already have gotten an approximation of for large q in (6.54). For small q, the 

volume of the momenta k, such that |k| < fcp < |k + q|, is proportional to q, as shown in 
Figure 6.3. Therefore, the integral in (6.53) transforms as follows 


—1 


E 


Q dk Aek.c 


|k|<A:p<|k+q| 


(27r)3 A4 q + z/2 


— i Ankp q d cos 6 cos 6 


; ,.2 q^/d + kpq^/d 

1 /vfT' o 


g^/2 + kpq cos 6 
{q‘^/2 + kpq cos 9)'^ + j^2 

q‘^ + 0 [q^) 


-ikp 


{q‘^j2) + jv2 


(6.56) 


where we used that Aek,q = (?^/2 + kpq cos 9 and that u is large compared to kpq- This is the 
same behavior of xo(Q')fo) as for large momenta q. For intermediate q, where the integration 
volume for k is roughly independent of q, the k integration merely averages the contributions 
to xo- For large u we retrieve 


—1 


E 


|k|<fcF<|k4-q| 


n dk Aek.q 

(27r)3 q + zy2 


j /I q"^/2 + kpq cos 9 

—-— / d cos 9 - K - 

3 J-i {q‘^/2 + kpq cos 9)^ + z^2 




.3 Q^ + OjQ) 

(g2/2)2 + i .2 


(6.57) 


Since the integrand has the same form for large ly in all cases, we can use it in the integral 
over the whole domain of q, getting 


Jo I ~ +Xo(g,fo)F(g)| 

4VF(z/2 + i)3/4-4i/2_3 _ 11 5 1 M 

16Z/5/2 ^ 192Z/9/2 + \^;,13/2 J ’ 

which is expanded in the variable u = 1 /j^ at n = 0. The other orders in the terms (6.56) and 
(6.57) yield a leading order term of the form 0(l/z/^) for the respective integration region 
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of q and can therefore be neglected for sufficiently large frequencies i/. These terms can, 
however, contribute to the next-to-leading order. Finally, we can insert this expansion in the 
imaginary frequency integration to estimate the missing RPA energy per electron AE^^^ 
when truncating the frequency integration at some finite but large 






+ 0 



3/2 

^^max 


+ o 



(6.58) 


Summary 


The Random Phase Approximation is the sum of all ring diagrams to infinite order. Mo¬ 
mentum conservation dictates that every Coulomb interaction in a ring diagram mediates 
the same momentum giving rise to a divergence in n-th order. This is the strongest 

divergence possible in n-th order rendering the ring diagrams the most important contribu¬ 
tion for low momenta, i.e. at long distances or in the high density regime. The divergence 
of each diagram when integrating over low momenta is referred to as infrared catastrophe. 
Evaluating the sum over all orders before integrating over the mediated momenta turns the 
Ijq^'^ divergence of each order into a log(l -|- I/q^) divergence, which yields a finite result in 
the subsequent momentum integration and solves the infrared catastrophe. 

Within the framework of many-body perturbation theory, discussed in Chapter 5, the 
Random Phase Approximation can be readily derived using the independent particle polariz¬ 
ability 0 = Xo as a building block. In the frequency domain this can be done using Feynman 
diagrams where the rotational symmetry of the ring diagrams gives rise to the factors in the 
expansion of the RPA energy 

The RPA can also be derived lii^the frequency domain within the Adiabatic Connection 
(AC) arriving at a formally equivalent result for very different reasons. In many-body per¬ 
turbation theory the perturbation is slowly introduced to the system leaving it in its ground 
state. The sum over all connected diagrams, respecting their symmetry, yields the total corre¬ 
lation energy. In the Adiabatic Connection the correlation energy is retrieved from averaging 
the potential energy over the coupling strength A: 


dA A 


({Lr-i{Lr 




In the AC the polarizability xa is the key quantity of interest rather than connected diagrams. 
Therefore, there are no symmetries to consider and connected (closed) diagrams should be 
avoided for depicting the RPA within the Adiabatic Connection. It is important to remark 
that the AC derivation is tailored to a DFT reference system assuming an exact density of 
the reference for the Adiabatic Connection and an exact density response for the Dyson-like 
equation of the polarizability. 

Finally, the RPA can also be derived in the time domain using Goldstone diagrams and 
the direct ring Coupled Cluster Doubles amplitudes 
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In each iteration of the amplitude equation ring diagrams are added in all possible ways in a 
time ordered manner. The simple fact that rings have a left and a right side when building 
them bottom to top requires 4 open connections of this building block according to the left 
and the right particle/hole pair. Although this approach is not as efficient as calculating the 
RPA in the imaginary frequency domain, it is easy to include a larger set of diagrams once 
the amplitudes are found. After all, Wick’s theorem requires all contractions to be considered 
not just the ring diagrams. The additional set of diagrams available given the direct ring 
Coupled Cluster Doubles amplitudes is called Second Order Screened Exchange diagrams. 
They represent the lowest order correction to the Random Phase Approximation. 




Chapter 7 

Second Order Screened Exchange 


The Random Phase Approximation (RPA) improves considerably on Hartree-Fock or Den¬ 
sity Functional Theory results. R is capable of describing van der Waals interactions and 
works well in a large variety of chemical environments. However, it is biased and tends to 
overestimate the negative correlation energy. In the previous chapter we have seen that for 
the uniform electron gas but this has also been shown for various solids and molecules. It is 
not surprising that the Random Phase Approximation shows an error. Wick’s theorem states 
that all contractions should be considered not just those forming the ring diagrams. That the 
RPA exhibits a systematic error, however, indicates a general reason behind this error which 
can serve as a guide to the next class of diagrams partially correcting RPA’s bias. 

In second order the overestimation of the correlation energy is more evident. According 
Figure 5.5 the second order (MP2) direct term is negative and reads 

'>(2-1-2) ^ Dj ab 

2 £i + Sj — £a — £b 
tjab 

where the Fermion sign is explicitly given depending on the number of loops and holes accord¬ 
ing to (5.38). This term contains contributions in violations of the Pauli exclusion principle, 
for instance as illustrated below on the left, where the state i occurs twice at the same in¬ 
stance in time. Such contributions are canceled exactly by the respective exchange diagram 
where the two offending states are crossed by anti-symmetrizing the affected interaction as 
shown on the right. The resulting diagram has one loop less giving an opposite sign; 



i)(2+2) 1 ^ _ /(_i)(i+2) 1 

2 £i £i — £a — £b i 2 Ej -|- ~ J 

As a consequence, violating contributions would be canceled if all contractions were consid¬ 
ered. Ignoring this exchange diagram leaves the violating negative contributions uncanceled 
resulting in a negative error. For the RPA the sign of the ring diagrams alternates with the 
order. However, the contributions of each order decay such that the error of the lowest order 
dominates, which is negative. 


51 






52 


CHAPTER 7. SECOND ORDER SCREENED EXCHANGE 



Figure 7.1: Exchanging only the last interaction in the ring diagrams leads to the Second 
Order Screened Exchange (SOSEX) correction to the Random Phase Approximation. 

The systematic error can be alleviated by including exchange diagrams where violations 
of the Pauli exclusion principle can occur. The lowest order correction to the Random Phase 
Approximation anti-symmetrizes only one interaction of the RPA ring diagrams. The correc¬ 
tion is termed Second Order Screened Exchange (SOSEX) if only the last interaction in time 
is anti-symmetrized. Eigure 7.1 shows the additional diagrams included in the SOSEX. 

7.1 SOSEX from Direct Ring Coupled Cluster Doubles 

The lowest order correction anti-symmetrizes only one interaction of the RPA ring diagrams. 
This could be any interaction and not necessarily just the last one. However, when calculating 
the Random Phase Approximation using the direct ring Coupled Cluster Doubles amplitudes 
tfj, as discussed in Section 6.2, calculating the SOSEX corrections comes with hardly any 
additional costs. We simply close the amplitudes once with a direct interaction and once 
with two indices swapped. Respecting the^Fermion sign the RPA-I-SOSEX energy is then 
given by 



The Second Order Screened Exchange correction to the RPA was introduced by (Ereeman 
1977) who applied it to the uniform electron gas. The term SOSEX was later coined by 
(Griineis et al. 2009) who studied this correction also for solids. Figure 7.2 shows the Second 
Order Screened Exchange correction per electron for the uniform electron gas as calculated by 
Freeman. The error of the Random Phase Approximation with respect to Quantum Monte- 
Carlo (QMC) results of (Ceperley and Alder 1980) fitted by (Perdew and Zunger 1981) is also 
given. Remarkably and certainly fortuitously RPA-I-SOSEX matches the QMC results at ~ 
5, which is in the density region of real metals. Anti-symmetrizing only the last interaction 
seems gives just the right correction to the RPA in the uniform electron gas. However, from a 
strictly ab-initio point of view, there is no reason to restrict the exchange diagrams to those 
where only the last interaction is anti-symmetrized, except technical convenience when having 
the drCCD amplitudes at hand. 

Anti-symmetrizing each but still only one interaction in the ring diagrams gives worse 
agreement with QMC in the high density regime r* < 8 but better agreement for low densities, 
where correlation effects are stronger. This is shown in Section 8.2. Note that in a spin- 
polarized system SOSEX is less fortunate. While still improving on RPA in the density range 
of interest, it cancels the RPA energy in the low density limit, as discussed in Subsection 
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Figure 7.2: The Second Order Screened Exchange energy per electron obtained by (Freeman 
1977) for the uniform electron gas. It is compared to the error of the Random Phase Approx¬ 
imation with respect to Quantum Monte-Carlo (QMC) calculations by (Ceperley and Alder 
1980) fitted by (Perdew and Zunger 1981). At « 5 RPA-I-SOSEX fortuitously matches 
QMC results. The Adiabatic Connection-SOSEX according to (7.16), which is discussed 
later, is also shown here. 


Ts 

[a.u.] 

(Ec-E^pA) 

[mEhN] 

^SOSEX 

[mEhN] 

^AC-SOSEX 

[mEh N] ± 

1 

19.167 

19.680 

19.832 

0.009 

2 

16.710 

17.560 

17.780 

0.003 

3 

15.545 

16.090 

16.342 

0.003 

4 

14.752 

14.970 

15.237 

0.003 

5 

14.131 

14.070 

14.343 

0.003 

6 

13.613 

13.320 

13.595 

0.003 

7 

13.165 

12.670 

12.955 

0.003 

8 

12.768 

12.120 

12.398 

0.003 

9 

12.413 

11.620 

11.906 

0.003 

10 

12.090 

11.190 

11.466 

0.003 

12 

11.521 

— 

10.712 

0.003 

15 

10.810 

— 

9.806 

0.003 

20 

9.884 

— 

8.679 

0.003 

30 

8.582 

— 

7.201 

0.003 

40 

7.685 

— 

6.246 

0.003 

50 

7.014 

— 

5.564 

0.003 


Table 7.1: Data of Figure 7.2 including low densities. 
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8.2.3. 

One can also include full anti-symmetrization in the amplitude equation, including terms 
such as 



which lead to the full Coupled Cluster Doubles (CCD) amplitude equation. However, calcu¬ 
lating the CCD amplitudes requires 0{N^) operations while the direct ring CCD amplitudes 
of RPA and SOSEX can be computed in 0{N^) steps (Scuseria and Schaefer 1989). A major 
drawback of either method is the large memory requirement scaling like 0{N'^) since the 
amplitudes need to be stored for iterating the amplitude equation. 


7.2 Adiabatic Connection-SOSEX 


In contrast to calculating the drCCD amplitudes, evaluating the Random Phase Approxima¬ 
tion in the frequency domain only requires 0{N‘^) memory. Calculating an exchange correc¬ 
tion based on two point quantities, such as the independent particle polarizability 0 = Xo 

would therefore be favorable. Angyan et al. suggested an approximation to the drCCD 
SOSEX within the framework of the Adiabatic Connection that can be implemented with a 
memory usage of 0{N‘^), which we will outline here. 

Within the Adiabatic Connection one can define a screened interaction similar to the one 
defined in many-body perturbation theory in (6.7) 



WxM = AV + A2VXo(iz^)V + A3VXo(izz)VXo(m)V + 


using the matrix notation introduced in Section 6.1. The only difference to (6.7) is the 
dependence on the coupling strength A. Note that we explicitly write the coupling strength 
in all diagrams within the Adiabatic Connection to make a clear distinction from the diagrams 
within many-body perturbation theory discussed in Chapter 5. 

We can now define the coupling strength averaged screened interaction 

W(m)= [ dAWA(m) = Jv +JvXo(ii/)V + ... 

Jo 2 6 


and write the RPA correlation energy found in (6.34) in terms of W: 


^RPA ^ / zAtr 

^ 2 / 27r 


1 fdiy 


{Xo(m)VXo(m)W(m)}. 


(7.3) 


Next, we insert the independent particle polarizability given in (6.28), 


V r ^ f>a(x)C(x')V’i(x')V'r(x) , V’i(x)V'r(x')V’a(x')V’a(x) 

Xoxx'fl^J ^ “ 2^ ' -^-^ 


Ea- Ei- IV 




lU 
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into (7.3), giving 

pRPA 1 f du 

- -o / ^ 2^ 




+ 




ijab 


(Sa - Ei- iu){eb - Ej - iu) (Ea - Ei - iu){Eb - Ej + ip) 


+ 


+ 


{Ea - Ei + il') {Eb - Ej - ip) {Ea - Ei + ip)(efe - Ej + il^) J ' 
where we write analogous to the matrix elements of the Coulomb operator; 

dxi dx2 V’p(xi)V'g(x2) WxiX2(i'^)'0r(x2)V's(xi) . 




(7.4) 


Although in the Adiabatic Connection the polarizability is the central quantity of interest 
rather than connected diagrams, we can use similar diagrams to depict the terms in (7.4): 



The diagrams are drawn such that the imaginary frequency p goes from right to left on the 
Coulomb interaction as indicated by the arrow. 

For real valued spin-orbitals V’p(x) the Coulomb integrals exhibit time reversal symmetry 
at each vertex such that . The same holds for the screened Coulomb 

integrals since the independent particle polarizability Xo(ip) is real valued. This simplifies 
(7.4) to 


pRPA ^ ^ 


dp 


E 

ijab 




IP 


)fA 


IP 


with /ja(ip) = 


2(£a Ei') 
{Ea - ej)2 -h p2 


In this form, the frequency dependent RPA energy expression bears resemblance to the drCCD 
RPA expression ^ define the AC-SOSEX by anti-symmetrizing the 

Coulomb interaction V in analogy to the drCCD SOSEX expression, arriving at 


^AC-SOSEX 



^ yaby^ij Mfia{y)fjb{y) ■ 


(7.6) 


Above equation still requires 0{N^) of memory from the two interactions and Wij. 
We can, however, transform it back into the position basis as discussed in Section 5.7, giving 

^AC-SOSEX , _1 y 

defining the imaginary frequency dependent exchange polarizability for the AC-SOSEX 
" “// dX3dx4 ^ ^V^*(x4)V^j(xi)V>;(xi)V^a(x3)/ia(ip) 

ia 

^V’j(x3)V'i(x2)V’fe(x2)V'b(x4)/i&(ip) • (7.7) 

jb 
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Note that the exchange polarizability is defined negative since it contains only one Fermion 
loop. We can give a closed form for the coupling strength averaged screened interaction 
finding a matrix function with the same Taylor expansion. This yields the final expression 
for the AC-SOSEX: 

^AC-SOSEX^^l J ^tr|pAC(ij,)(^Xo(m)VXo(m))“'(log(l-Xo(m)v)+Xo(m)v)| . 

(7.8) 

The exchange polarizability is a quantity depending on two positions. Thus, evaluating 
the AC-SOSEX as given above requires only 0{N‘^) of memory rather than 0{N^), greatly 
broadening the applicability of the AC-SOSEX to larger systems. However, calculating the 
exchange polarizability still requires 0{N^) steps, which is equally time consuming as calcu¬ 
lating the SOSEX from the direct ring Coupled Cluster Doubles amplitudes. In the Random 
Phase Approximation the respective polarizability analogous to P^'^(iz/) simply factors into 
Xo(R)VXo (ii/), allowing for an evaluation in only 0{N^) steps. This can not be done for the 
AC-SOSEX since X 3 and X 4 occur in both sums in (7.7). Reducing the memory consumption 
from 0{N‘^) to 0{N‘^) is still an important improvement since it is easier to allocate more 
CPUs to a calculation than it is to allocate more memory per CPU. 

7.3 Difference between drCCD SOSEX and AC-SOSEX 

Despite the similarity of (7.6) and the expression for the drCCD SOSEX (7.1) they are 
only identical in second order but not beyond. First, the drCCD amplitudes V.V are 
constructed monotonous in time according to (6.13)-(6.16), which guarantees that the closing 
Coulomb interaction in 



is indeed the last interaction in time. In contrast, the averaged screened interaction W(iz/) 
contains interactions that reach both, in the past and in the future. Thus, the left diagram 
shown in Figure 7.3 is contained in the AC-SOSEX while it is not contained in the drCCD 
SOSEX. Furthermore, anti-symmetrizing the first and the last case in (7.5) yields terms that 
have no correspondence in many-body perturbation theory. The AC-SOSEX introduces anti- 
symmetrization by simply swapping i and j at the unscreened interaction V. For the last 
term this yields for example 

vSwi. 

This term, however, contradicts the requirement of many-body perturbation theory that 
upper indices can only match lower indices and vice-versa since upper and lower indices 
refer to creation and annihilation operators, respectively. This term can therefore not be 
drawn diagrammatically in the usual manner such that the upper/lower indices correspond to 
outgoing/incoming connections. We can, however, draw the term respecting the propagation 
direction of particles and holes. The resulting diagram is depicted on the right of Figure 7.3. 
It exhibits a particle a turning into a hole j at the left vertex of the Coulomb interaction 
and a hole i turning into a particle h at the right vertex. This diagram will be termed 
swapped ladder diagram since it resembles a particle-hole ladder diagram where b and j are 
swapped. Employing the same notion of diagrams as in (7.5), the AC-SOSEX can be depicted 
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Figure 7.3: Diagrams contained in the AC-SOSEX that are not part of the SOSEX based on 
the drCCD amplitudes. The factors from the coupling strength integration are omitted. The 
right diagram shows a particle turning into a hole and vice-versa, which has no correspondence 
in many-body perturbation theory. 


diagrammatically by 



7.3.1 Diamond C(A4) 

We can study the swapped ladder diagram for a small test system consisting of diamond C(A4) 
with 128 states. It is based on a Hartree-Eock reference with 2x2x2 fc-points in a primitive 
cell comprising 2 atoms with 4 electrons per atom. Although the system is only coarsely 
described by 128 states it serves well as a benchmark for individual diagrams beyond second 
order. In this finite band gap system we can limit the order of the AC-SOSEX diagrams 
from (7.9) solely to third order while still getting finite results. This excludes all diagrams 
that are not contained in the drCCD SOSEX except the third order swapped ladder diagram. 
The AC-SOSEX expression from (7.9) expands in third order to all possible permutations of 
the interaction times, analogous to Eigure 5.9. Including the factor 1/3 from the coupling 
strength integration this gives 



(7.10) 


The lower diagrams are equivalent to the upper diagrams since each can be continuously 
deformed into the respective upper diagram without changing the order of the Coulomb 
interactions. The left two and the right two diagrams are identical due to time reversal 
symmetry. This leaves two diagrams with distinct energies to evaluate, the exchange diagram 
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and the swapped ladder diagram; 




/_]^'| 2+3 _ ''ab'^ik '^jc _ 

i^abc A Sk — Sa — + £j ~ £a ~ £fe) 


—1.151 mEh N 


(7.11) 


(- 1 ) 


2+3 


ai ik be 


^t^bc {A + eu-ea-ec){e^ 


+ £k — £b — 


—1.141 uiEh N 


(7.12) 


Although the denominators differ, they yield almost the same energy per electron. Thus, the 
AC-SOSEX energy for the test system in third order hardly differs from the drCCD SOSEX 
energy in third order. The latter is simply given by the exchange diagram in (7.11): 


£.AC-SOSEx( 3) = _ j*L(2 • 1.151 + 1.141) = -1.148 mE;,iV 

3 

EfOSEX(3) ^ ^ ^ 

In fourth order, the AC-SOSEX expression from (7.9) expands to 24 possible permutations 
of the interaction times, 12 of which are distinct Goldstone diagrams. This cancels the factor 
1/2 of (7.9). The factor from the coupling strength integration is 1/4, giving 



(7.13) 


where t.r. denotes the time reversed variants of the shown 6 diagrams. Note that the diagrams 
are drawn such that the permutations of the interaction times are evident and not according to 
minimal self intersection. The diagrams can be evaluated by iterating the doubles amplitudes 
employing only a subset of the direct ring Coupled Cluster Doubles amplitude equations. 
The following steps are for instance used to calculate the second diagram in (7.13) indicating 
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the employed part of the drCCD amplitude equation above each equals sign; 


v..y 

t(0) 


.ab(0) ( 6 ^ 3 ) 

'ij 


V 


Ei -\- Ej — - £b 




Evaluating the second row of diagrams in (7.13) requires two additional parts of the doubles 
amplitude equation that are not part of the drCCD amplitude equation: 


t = 0 


a i b J 

v__y 


t = 0 


a i b J 



(7.14) 


(7.15) 


In the case of the swapped ladder amplitude equation (7.14) there is one additional hole k but 
no additional loop, resulting in a negative sign. In the exchange amplitude equation (7.15) 
there are two additional holes and one additional loop also giving a negative sign. The latter 
is part of the full Coupled Cluster Doubles (CCD) amplitude equations while the swapped 
ladder equation only occurs in the AC-SOSEX. Table 7.2 lists the energies per electron of all 
diagrams of the AC-SOSEX in third and fourth order modulo time reversal symmetry. The 
AC-SOSEX energy in fourth order is thus 


_g.AC-SOSEX(4) 


1 

2 


(0.212 0.105 0.208 0.208 0.103 0.209) 


+Q.I)22uiEhN. 


The drCCD SOSEX energy is formed by the diagrams of the first row in Table 7.2 only, 
yielding a very similar result; 


EfOSExD) ^ Q 212 0.105 0.208 = -hO.525 mE,* N . 


The top row in Table 7.2 shows diagrams where the last interaction is anti-symmetrized 
while in the second row it is the second last interaction that is anti-symmetrized. The table 
indicates that the energy of the SOSEX or AC-SOSEX diagrams hardly depends on which of 
the interactions is anti-symmetrized when using the swapped ladder diagram instead of the 
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Table 7.2; Goldstone diagrams forming the third and the fourth order of the drCCD SOSEX 
and the AC-SOSEX correction. The drCCD SOSEX consists of diagrams of the first row 
only. A weighted average of the first and second row yields the AC-SOSEX energy. 


proper ladder diagram. Finally, the factor 1/n from the coupling strength integration simply 
averages the energies from anti-symmetrizing each of the n interactions, as shown in (7.10) for 
third order and in (7.13) for fourth order. Since the energies of the diagrams are very similar 
anti-symmetrizing each of the n interactions, the average is also very similar to the energy 
of the diagrams where only the last interaction is anti-symmetrized. The averaged energy 
forms the AC-SOSEX energy while the latter forms the drCCD SOSEX originally proposed 
by Freeman. This argument holds at least up to fourth order in the diamond test system and 
the truncation at fourth order describes the drCCD SOSEX already to an accuracy of 5%. 

7.3.2 Uniform electron gas 

In a metallic system it is not possible to truncate the AC-SOSEX expression at any finite order 
beyond the second order since all of the described higher order exchange diagrams diverge. It 
is, however, possible to study the AC-SOSEX diagrams of (7.9) numerically for the uniform 
electron gas (UEG) as the prototypical metal. In the UEG we let spin and momenta denote 
the states i,j, a and b occurring in the diagrams. For the third and the fourth diagram of 
(7.9) we choose for instance the following definition of momenta: 
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ki and —k 2 are required to be hole states while ki + q and —k 2 — q must be particle states. 
Thus, kj must be below and kj + q must be above the Fermi momentum kp for i G {1,2}. 
Note that in the swapped ladder diagram on the right the momenta —k 2 and —k 2 — q are 
dehned opposite to the propagation directions of the respective states, as indicated by the 
arrows next to the labels. The spins of all four states must be the same. In the UEG the sum 
over the states is replaced by the sum over the spin and integrals over all internal momenta 
according to (6.38), arriving at 


E. 


AC-SOSEX 


n 1 

~^N2 




dq 

{27ry 


E 


n^dkidka 


1 


(Aeki,q - ii^)(Ae-k2-q “ W 

1 


|ki|<fcF<|ki+q| 


+ 


(27r)f 


V (ki + k 2 + q)W (q, v) 


1 


(A^ki.q - m)(Ae_k2-q + 
1 


(Aeki,q + ij^)('^e-k2 -q “ (^eki,q + m)(A£_k2,-q + ii^) J 

with i G {1,2} and the single particle excitation energy Aek.,q = (k* + q)^/2 — k?/2. The 
same expression can be derived using imaginary frequency propagators defined in (6.44) and 
integrating out the two additional imaginary frequencies analogous to the derivation of xo in 
(6.45). 

Integrating out ki and k 2 turns out to be a tedious task for the above equation. Although 
there are closed expressions for the second and the third case they are overly complicated. For 
the swapped ladder diagrams in the first and the fourth case no such expressions were found. 
However, a straight-forward Monte-Carlo integration of ki and k 2 has proven to be sufficiently 
accurate when sampling the momenta kj with a probability density function (PDF) given by 

1 


PDF(kj 


oc 


Aeki,q ± i 




0 


for |kj| < fcp < |ki -|- q| 
otherwise. 


Figure 7.2 and Table 7.1 show the resulting AC-SOSEX energies as a function of density given 
by the Wigner-Seitz radius, r*. The uncertainties from the integrations are indicated by the 
error bars. Eor the Monte-Carlo integration of ki and k 2 a precision of 5 significant digits 
can be achieved with less than 30000 samples for each q and v, depending on momentum, 
frequency and density. The error from the momentum and imaginary frequency integration 
is of similar magnitude. 

The differences between the two SOSEX variants are not as small as for an isolating 
system but they are still below 3% for the density range with < 10. Eor lower densities no 
drCCD SOSEX reference values were found. 


Summary 


The Random Phase Approximation (RPA) systematically overestimates the negative correla¬ 
tion energy. This originates, at least partially, from violations of the Pauli exclusion principle 
in the ring diagrams of the RPA such as 
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By the merit of Wick’s theorem, violations of the Pauli exclusion principle do not have 
to be considered as long as all contractions of the occurring operators are included. The 
contractions correcting these violations are those from the respective exchange diagrams where 
the offending states are crossed by anti-symmetrizing an affected Coulomb interaction. 

Thus, the lowest order correction to the Random Phase Approximation anti-symmetrizes 
one Coulomb interaction occurring in the ring diagrams. If this interaction is the last inter¬ 
action in time the respective correction is termed Second Order Screened Exchange (SOSEX) 
containing the following diagrams: 



The SOSEX can only be computed from the direct ring Coupled Cluster Doubles (drCCD) 
amplitudes v.y since monotonicity in time is required to be able to anti-symmetrize only 
the last interaction. When calculating the RPA using the drCCD amplitudes the SOSEX 
can be easily computed with hardly any additional costs. However, calculating the drCCD 
amplitudes requires 0{N^) of memory limiting the applicability of RPA-I-SOSEX to rather 
small systems. 

To overcome the limitations regarding memory consumption the Adiabatic-Connection 
(AC) SOSEX can be used requiring only 0{N‘^) of memory. It yields very similar results 
compared to the SOSEX although it contains terms that have no correspondence in many- 
body perturbation theory where particles turn into holes and vice-versa herein called swapped 
ladder diagrams. This is due to the numerical oddity that it hardly matters which of the n 
occurring interactions in n-th order are anti-symmetrized as long as the swapped ladder 
diagrams are used rather than actual ladder diagrams. In third order this means for instance: 



such that 



Thus, the AC-SOSEX can be considered a recipe for imitating the SOSEX energy while 
reducing the memory requirements to O(iV^). 

Apart from technical convenience when having the drCCD amplitudes at hand, there 
is however no reason to limit the considered exchange diagrams to those where only the 
last interaction is anti-symmetrized. From an ab-initio point of view one should consider all 
diagrams where one interaction is anti-symmetrized and where violations of the Pauli principle 
can occur as the lowest order correction to the RPA. This is discussed in the next chapter. 



Chapter 8 


Adjacent Pairs Exchange 


The Second Order Screened Exchange (SOSEX) correction to the Random Phase Approxi¬ 
mation (RPA) arises from anti-symmetrizing only the last Coulomb interaction of RPA’s ring 
diagrams as shown in Eigure 7.1. When calculating the RPA using the direct ring Coupled 
Cluster Doubles (drCCD) amplitudes as introduced by (Ereeman 1977), this comes at 
virtually no extra costs and represents a natural choice for the lowest order correction to the 
RPA. However, when calculating the RPA in the frequency domain, which is more efficient, 
there is a priori no reason to choose this particular class of diagrams as the lowest order cor¬ 
rection to the RPA. Purthermore, it is not possible to evaluate the SOSEX diagrams directly 
since the last interaction in time cannot be explicitly addressed in the frequency domain. 
The AC-SOSEX approach, discussed in Section 7.2, offers an approximation but it contains 
swapped ladder diagrams that are not part of the many-body perturbation expansion. So 
the question remains, which diagrams, that are actually part of the many-body perturbation 
expansion, can be efficiently evaluated in the frequency domain and offer a good lowest order 
correction to the systematic error of the Random Phase Approximation. 


8.1 Drivation of the Adjacent Pairs Exchange correction 

As discussed in the beginning of Chapter 7, violations of the Pauli exclusion principle in RPA’s 
ring diagrams suggest to anti-symmetrize the Coulomb interactions wherever such a violation 
can occur. In lowest order it should be only one but not necessarily just the last Coulomb 
interaction to be anti-symmetrized. This is done by cutting out one Coulomb interaction 
including the adjacent pair bubbles from the RPA ring diagrams and anti-symmetrizing this 
interaction if that can correct for a violation of the Pauli exclusion principle. 


8.1.1 Two Sided Adjacent Pairs Exchange 

Two adjacent pair bubbles have four possible time orders with respect to the Coulomb in¬ 
teraction between them, shown in Eigure 8.1(a). In the first and in the last case an anti- 
symmetrization of the contained Coulomb interaction only cancels contributions where the 
same states i and a occur in consecutive bubbles, as illustrated in Eigure 8.1(b). In general, 
such contributions do not violate the Pauli exclusion principle. Pollowing this argument, we 
exclude these cases and study a correction to the Random Phase Approximation where both 
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Figure 8.1: (a) The possible time orders of a Coulomb interaction and its two adjacent pair 
bubbles, (b) Anti-symmetrizing the first and the last case only cancels contributions i and 
a as indicated here for example. In general, these contributions do not violate the Pauli 
exclusion principle. 



Vs [a.u.] 

Figure 8.2: The Two sided Adjacent Pairs Exchange (2APX) energy per electron for the 
uniform electron gas compared to the error of the Random Phase Approximation with respect 
to Quantum Monte-Carlo (QMC) calculations by (Ceperley and Alder 1980) fitted by (Perdew 
and Zunger 1981). Anti-symmetrizing more than one interaction in the RPA ring diagrams by 
the application of (8.1) worsens the accuracy with respect to Quantum Monte Carlo results. 
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remaining cases are anti-symmetrized to 



( 8 . 1 ) 


and inserted into the RPA ring diagrams. Since these diagrams are polarization parts with 
two open vertices the memory requirement for calculating this correction scales like 0{N‘^) 
with the system size N. This is equivalent to the memory requirement of the AC-SOSEX 
and considerably less than 0{N‘^) of the conventional SOSEX, calculated from the direct ring 
Coupled Cluster Doubles (drCCD) amplitudes. We term this correction two sided Adjacent 
Pairs Exchange (2APX) and its expansion in terms of Goldstone diagrams is given by 




+ 


t.r. + ... 



+ ... + 


( 8 . 2 ) 


The second row contains only the right contribution of (8.1) and t.r. refers to diagrams 
emerging from the second row by time reversal, containing only the left contribution of (8.1). 

The two sided Adjacent Pair Exchange correction differs from the Second Order Screened 
Exchange correction already in third order. Time reversal of any SOSEX diagram beyond 
second order gives a distinct diagram and that is contained in the two sided APX. Since third 
order is the lowest order where the two sided APX and SOSEX differ and since this order 
is in general negative, the two sided APX correction is expected to be less than the SOSEX 
correction. Evaluating the two sided APX for the Uniform Electron Gas, as discussed in 
Section 8.2, confirms this expectation. As shown in Eigure 8.2 the two sided APX considerably 
underestimates the desired energy correction, given by the difference of the RPA and Quantum 
Monte Carlo results. At low densities, where > 10, it even becomes negative, actually 
worsening the systematic error of RPA. 


The two sided Adjacent Pairs Exchange correction seems the most plausible lowest order 
correction to the Random Phase Approximation. However, despite improving on the RPA in 
the UEG for densities with r* < 10, the two sided APX does not offer a balanced correction 
to the RPA since it always underestimates, but never overestimates the missing correlation 
energy. One could include diagrams, where two, three or more of the RPA’s Coulomb in¬ 
teractions are anti-symmetrized, in the fashion discussed above, as the next orders of the 
correction. These corrections are still forming a ring and can thus be efficiently evaluated in 
the frequency domain, once the two sided Adjacent Pairs Exchange polarization part 
has been calculated. In terms of Feynman diagrams and propagator matrices, the corrections 
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Figure 8.3; Exchange of pair bubbles that are non-adjacent 


with one or two Coulomb interactions anti-symmetrized are given by 


^2APX _ 



^^2APX2) ^ 


with 



and where the superscript of the trace, the imaginary time integration - including the 

sign - and the imaginary time arguments have been omitted for brevity. The RPA screened 
interaction is given by (6.7). All diagrams posses reflection symmetry but note that 

the symmetry factor differs in the case where only one occurrence of is inserted into the 

ring diagrams compared to other case. This is due to a reflection symmetry introduced when 
closing with only one Coulomb interaction V. Figure 8.2 shows E^^^^ and 

according to (8.3) and (8.4). More insertions of P^^^^ into the ring diagrams of the Random 
Phase Approximation actually worsen the a accuracy of the two sided APX with respect to 
Quantum Monte Carlo results. More than two insertions of P^^^^ into the ring diagrams 
of the RPA offer no improvement either so none of the two sided APX approximations ever 
overestimates the missing correlation energy and thus none is balanced. 

We could investigate more complex exchange processes, for instance those correcting vio¬ 
lations of the Pauli exclusion principle of pair bubbles in the RPA that are not adjacent, as 
sketched in Figure 8.3. However, the exclusion principle is only violated if the two effected 
pair bubbles propagate at overlapping times. This would require a time order which can 
only be provided when evaluating the Random Phase Approximation from the direct ring 
Coupled Cluster Doubles (drCCD) amplitudes, loosing the advantage of the reduced memory 
requirements of an RPA implementation in the frequency domain. 


\ + P-VXoV + PxVXoVXoV + .. .^ (8.3) 


^ (PxV + PxVXoV + P.VXoVXoV + ...)" (8.4) 
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Figure 8.4: The two sided Adjacent Pairs Exchange correction introduces new violations of 
the Pauli exclusion principle, as shown here in third order in the case of two identical hole 
indices. One of the two exchange diagrams exactly cancels the offending contributions of the 
RPA diagram, indicated by parenthesis and blue index labels. The other exchange diagram, 
however, introduces new violations shown by red index labels. 


8.1.2 Adjacent Pairs Exchange 

The unbalanced performance of the two sided APX and its higher order variants rises the 
question whether really all Coulomb interactions should be anti-symmetrized if that can 
correct for violations of the Pauli exclusion principle. It turns out that the two sided APX, 
while indeed correcting for all violations occurring in adjacent pairs, introduces new violations. 
In third order this is most apparent and illustrated in Figure 8.4. The lower/upper row shows 
the case where the lower/upper two pair bubbles of the RPA propagate in the same hole state 
i. The diagram shown with blue index labels exchanges these two propagators and exactly 
cancels the offending contributions of the RPA. This is indicated by parenthesis around the 
RPA and the respective exchange diagram. The other diagram, however, introduces new 
violations to the Pauli exclusion principle, shown by the red index labels. 

The key issue in this case is that the leftmost and longest pair bubble of the RPA diagram 
is exchanged in both diagrams of the two sided APX corrections. While this guarantees that 
all possible violating contributions are canceled, it always introduces new violations. For a 
more balanced correction we therefore require that each pair bubble is exchanged at most 
once. The simplest correction satisfying this requirement is the (single sided) Adjacent Pairs 
Exchange (APX) correction, where only one case of (8.1) is contained. For a system with 
time reflection symmetry it is irrelevant which of the cases we choose and without loss of 
generality we choose to anti-symmetrize adjacent pair bubbles in the third case of time orders 
shown in Figure 8.1(a): 



(8.5) 


In terms of Feynman diagrams and the matrix notation of the propagators, introduced in 
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Section 6.1, the APX correction is then given by 




{pAP^w} 


} 


( 8 . 6 ) 


with 




and where the imaginary time arguments as well as the superscript of Px in the diagrams 
have been omitted for brevity. All diagrams exhibit a single reflection symmetry but note 
that, unlike in the two sided APX case, there is no additional symmetry introduced when 
closing Px with one Coulomb interaction V since the Px^^ contains only one of the two 
time orders contained in Px"^P^. Despite the ring form of the APX diagrams beyond second 
order, none of them has a rotational symmetry in contrast to the respective RPA diagrams. 
This simplifies the sum over all orders of the perturbation compared to the RPA since all 
orders have the same factor. We can use the infinite sum of a geometric series to give an 
explicit form for the APX energy: 

= ^tr|pAPX(i^)v(l-Xo(R)v)''| , (8.7) 


where the imaginary frequency dependent exchange polarizability for the APX contains only 
one of the four time orders contained in Px*^, which was given in (7.7): 

= - [[ dx3dx4 ^ ^V^*(x4)V^»(xi)V>;(xi)V:a(x3) ■ 

J J ^-3 ^4 “T 

la 

^V’j(^3)V'i(x2)V’fe(x2)V'b(x4)- _ ^ . • (8.8) 


8.2 APX for the uniform electron gas 

Evaluating the Adjacent Pairs Exchange (APX) energy for the Uniform Electron Gas (UEG) 
is very similar to evaluating the AC-SOSEX energy according to (7.16). In contrast to the AC- 
SOSEX expression, only one of the four time orders of two adjacent pair bubbles are contained 
in the APX. In the chosen order, the Goulomb interaction of the exchange polarizability Px 
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occurs after both open vertices of 

ki + k2 + q 



Furthermore, the APX uses the screened interaction W, represented by the double wiggly 
line, as given by (6.7) rather than the coupling strength averaged screened interaction W, 
since all orders of the expansion have the same factor. The APX correction to the Random 
Phase Approximation per electron is thus given by 


E. 


■APX 


n 1 r dv 
~ '^N2 / ^ 


dq 


( 2 vr; 


.3 E 


|ki|<fcF<|ki+q| 


f^Mkidka 

(27r)® 


V (ki + k 2 + q}W (q, u) 


(^^ki,q T 1^)k2,—q 




(8.9) 


with i G {1, 2} and the single particle excitation energy Aeki,q = (1^* + t)^/2 —kf/2 . As in the 
case of the AC-SOSEX, the above expression can be evaluated by a Monte-Carlo integration 
for ki and k 2 using a probability density function (PDF) given by 


PDF(ki) oc 


1 


Aek„q ± ii' 


0 


for |kj| < fep < |ki +q| , 
otherwise. 


The momentum q and the imaginary frequency p can be integrated using a Gauss-Kronrod 
rule, analogous to the evaluation of the Random Phase Approximation. However, the asymp¬ 
totic behavior of the APX energy for large imaginary frequencies differs from that of the RPA 
and the AC-SOSEX, such that a different variable transform for integrating the frequency tail 
must be used. The asymptotic behavior for large frequencies is discussed Subsection 8.2.2. 

Eor the lowest order of the APX, the accuracy of the numerical integrations can be bench- 
marked against the MP2 exchange energy, which is independent of the density and it is known 
analytically from the work of (Onsager, Mittag, and Stephen 1966): 


log(2) 3 C(3) 


(2vr) 


0.02417915891814441Eft N. 


For the Monte-Carlo integration of ki and k 2 a precision of 5 significant digits can be achieved 
with less than 30000 samples for each q and v, depending on momentum, frequency and 
density. The error from the momentum and imaginary frequency integration is of similar 
magnitude. Figure 8.5 and Table 8.1 show the resulting APX energies as a function of density 
given by the Wigner-Seitz radius, r^. The uncertainties from the integrations are indicated 
by the error bars. 

^In the given diagram, the time order of the vertices is only relevant with respect to the Coulomb interaction. 
It is therefore neither a Goldstone nor a Feynman diagram. The proper Feynman diagram of APX is the 
leftmost diagram of (8.6). 



70 


CHAPTER 8. ADJACENT PAIRS EXCHANGE 



Figure 8.5: The Adjacent Pairs Exchange (APX) energy per electron for the uniform electron 
gas compared to the error of the Random Phase Approximation with respect to Quantum 
Monte-Carlo (QMC) calculations by (Ceperley and Alder 1980) fitted by (Perdew and Zunger 
1981). RPA-I-APX also fortuitously matches the QMC results but only at ~ 10. However, 
in the low density regime, where correlation is stronger, APX is considerably closer to the 
QMC results than AC-SOSEX. For results of the spin-polarized uniform electron gas see 
Subsection 8.2.3. 


Ts 

[a.u.] 

(Ec - Ef^^) 
[uiEhN] 

[mE,, N] ± 

^AC-SOSEX 

[mE^AT] ± 

1 

19.167 

19.956 

0.005 

19.832 

0.009 

2 

16.710 

18.012 

0.005 

17.780 

0.003 

3 

15.545 

16.672 

0.005 

16.342 

0.003 

4 

14.752 

15.649 

0.005 

15.237 

0.003 

5 

14.131 

14.825 

0.004 

14.343 

0.003 

6 

13.613 

14.139 

0.004 

13.595 

0.003 

7 

13.165 

13.553 

0.004 

12.955 

0.003 

8 

12.768 

13.043 

0.004 

12.398 

0.003 

9 

12.413 

12.594 

0.004 

11.906 

0.003 

10 

12.090 

12.194 

0.004 

11.466 

0.003 

12 

11.521 

11.506 

0.004 

10.712 

0.003 

15 

10.810 

10.680 

0.004 

9.806 

0.003 

20 

9.884 

9.650 

0.004 

8.679 

0.003 

30 

8.582 

8.289 

0.004 

7.201 

0.003 

40 

7.685 

7.400 

0.004 

6.246 

0.003 

50 

7.014 

6.757 

0.004 

5.564 

0.003 


Table 8.1: Data of Figure 8.5 including low densities. 
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8.2.1 Large momentum behavior 

In a solid or in a molecule, the finite resolution of the DFT or the Hartree-Fock reference 
imposes an upper limit Gmax on up to where the exchange polarizability Px and the inde¬ 
pendent particle polarizability Xq can be evaluated. Assuming that the system is sufficiently 
homogeneous at that resolution the asymptotic behavior for large momenta in the Uniform 
Electron Gas can be used to extrapolate results retrieved at finite Gmax to the limit of an 
infinite basis set, where Gmax —> oo. 

We can write the APX correction in the UEG from (8.9) in terms an exchange polariz¬ 
ability analogous to the general APX expression (8.7); 


;nAPX _ ___ 

N2 


= -Yl 



|ki|<fcp<|ki4-q| 

( 8 . 11 ) 

To get an approximation P^^^ for large momenta q we can trivially integrate ki and k 2 , 
since |kj| < /cp ( 7 , getting 




1 


1 


(g 2 / 2 -|-iz/)(g 2/2 — iz/) 


( 8 . 12 ) 


Inserting this and the approximation of Xo{<lA^)V{q) for large q, given in (6.54), into the 
APX energy expression, we can integrate p for a given, large momentum q: 

{q'^l2f + v‘^ ) - 4 

The missing energy of the APX correction when truncating the momentum integra¬ 

tion at a finite but large momentum Gmax is thus 


di/ 1 1 

27r q'^ (( 7^/2 -I- \v){q^ 12 — izz) 


1 - 


AE, 


APX 


(G„ 


'G„ 


- 4 


as 


G3, 


-I- a^ 


GL 


+ 0 


Gii 


(8.13) 


where we have expanded the integrand in the variable u = 1/q at u = 0. Hence, the APX 
energy and, similarly, the AC-SOSEX energy have the same asymptotic behavior with respect 
to large momenta q as the Random Phase Approximation. However, the convergence with 
respect to Gmax is slower since the exchange corrections are more short ranged, which may 
require the use of more than the leading order term of the Taylor expansion for an accurate 
extrapolation to the infinite basis set limit. 


8.2.2 Large imaginary frequency behavior 

For the frequency integration of metallic systems, such as the Uniform Electron Gas, knowl¬ 
edge about the asymptotic behavior of the APX correction with respect to large imaginary 
frequencies is important for choosing an appropriate variable transform for the tail. 

Eor large momenta q, we can use the approximation of P^^^ in (8.12), gotten in the 
previous subsection. Considering large frequencies we can still use the same approximation 
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Figure 8.6: Cross section of the set of momenta kj such that |kj| < kp < |kj + q| for excitation 
momenta q < 2. The momentum of the Coulomb interaction ki+k 2 +q vanishes for the shown 
ki and k 2 . The vicinity of the singular ki is magnified on the left, showing the employed 
coordinates for the integration of this regions and a volume element in form of a ring with 
area 27r|q|ui. 


for intermediate q > 2/cf since the denominator of the propagator V (ki + k 2 + q) is non¬ 
vanishing and the volume of the integration region of ki and k 2 is independent of q. In this 
case the integration of ki and k 2 merely averages the contributions. For q < 2/cf, the volume 
of the integration region of ki and k 2 depends on q and there are also contributions where 
the propagator C(ki -|- k 2 -|- q) becomes singular, as shown in Figure 8.6. In this case, we 
split the integration into two parts. In the first part both, ki and k 2 are in the vicinity of 
the singular contribution, shown for ki on the left of Figure 8.6. In the remaining part either 
ki, k 2 or both are away from the singular contribution. For the first part, we can transform 
the integration of ki and k 2 into an integration of the respective distances ui and U 2 from 
the singular contribution in the direction of q. In the remaining part, we approximate the 
integral assuming that ki -|- k 2 average out and taking into consideration that the integration 
volume of each kj scales like q^ + q'^. We get 


|ki|<fcF<|ki+q| 


IlMkidk2 

(27r)® 


V (ki -I- k 2 -I- q) 


k-p—ql2 


duidu] 


(27r)^g^uiU2 
{ui + U2Y 


+ 


{a^q-^ + 02 ^) 


2\2 


<f + 0 {q^) . 


(8.14) 


Since q < 2 and u is large, we can ignore the propagators of P^^^ containing the imaginary 
frequency for the integration of ki and k 2 . 

We can now insert the respective approximation of P^^^ for small q, intermediate q 
and large q into the expression for the APX energy (8.10) and integrate q. For xo we use 
the approximations (6.56), (6.57) and (6.54) for small, intermediate and large momenta q, 
respectively. The results are long and apart from their expansion in 12 of no particular interest. 
For small and intermediate q the integration of q yields A O and for large q we 

get l/i/^-|-O (l/z/^/^). Finally, we can insert this expansion in the imaginary frequency 
integration of the APX energy and estimate the missing energy per electron A£'^^^(i/max) 



8.2. APX FOR THE UNIFORM ELECTRON GAS 


73 



Figure 8.7: Numerical comparison of the missing correlation energy per electron for the 
Random Phase Approximation and its corrections when truncating the imaginary frequency 
integration at I'max- t'max is given in units of fcp and kp = 1 a.u. 

when truncating the integration at some finite but large imaginary frequency 



The asymptotic behavior of the APX energy differs from that of the Random Phase Approxi¬ 
mation and from that of the AC-SOSEX, which can be derived in an analogous fashion. This 
originates from the imaginary frequency behavior of containing only one of the four 

time orders contained in the RPA and in the AC-SOSEX. Figure 8.7 compares the asymptotic 
behavior of RPA, AC-SOSEX and APX for the Uniform Electron Gas numerically. 

8.2.3 Spin-polarized Uniform Electron Gas 

We can readily evaluate the RPA and the Adjacent Pairs Exchange correction for the spin- 
polarized Uniform Electron Gas using only one spin in the sum over all spins occurring 
in the expression of ii^) in (8.11) and Xo(q)i^) ™ (6.46). Considering that kp also 

depends on wheter the UEG is spin-polarized or not according to (6.39), yields correlation 
energies of RPA-I-APX shown in Figure 8.8 and Table 8.2. For comparison, the correlation 
energy of RPA-|-AC-SOSEX according to (7.16) is also given. 

Unlike in the non-spin-polarized case, neither of the SOSEX variants offers a balanced 
correction to the Random Phase Approximation, overestimating the missing energy for the 
entire range of densities. The accuracy with respect to Quantum Monte Carlo results also 
worsens for low densities. For low densities we can assume ki-|-k 2 -|-q ~ q since |kj| < kp <C 1. 
In this limit, two adjacent RPA bubbles only differ from the respective exchange diagram in 
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Figure 8.8: The Adjacent Pairs Exchange (APX) energy per electron for the spin-polarized 
uniform electron gas compared to the error of the Random Phase Approximation with respect 
to Quantum Monte-Carlo (QMC) calculations by (Ceperley and Alder 1980) fitted by (Perdew 
and Zunger 1981). For low densities, anti-symmetrization of a ring diagram merely changes its 
sign, such that SOSEX or its variants simply remove the correlation energy already captured 
by the RPA. 


Ts 

[a.u.] 

(Ec - 
[uiEhN] 

EAPX 

[mEh A^] ± 

^AC-SOSEX 

[mEh A^] ± 

1 

20.192 

21.809 

0.005 

21.746 

0.033 

2 

18.326 

20.497 

0.005 

20.394 

0.017 

3 

17.131 

19.511 

0.005 

19.374 

0.006 

4 

16.218 

18.712 

0.005 

18.525 

0.003 

5 

15.472 

18.037 

0.005 

17.805 

0.003 

6 

14.840 

17.452 

0.005 

17.179 

0.003 

7 

14.293 

16.936 

0.005 

16.627 

0.003 

8 

13.809 

16.475 

0.005 

16.130 

0.003 

9 

13.377 

16.058 

0.005 

15.680 

0.003 

10 

12.987 

15.678 

0.005 

15.268 

0.003 

12 

12.307 

15.007 

0.005 

14.541 

0.003 

15 

11.469 

14.169 

0.004 

13.628 

0.003 

20 

10.398 

13.073 

0.004 

12.431 

0.003 

30 

8.932 

11.536 

0.004 

10.745 

0.003 

40 

7.943 

10.474 

0.004 

9.580 

0.003 

50 

7.215 

9.678 

0.004 

8.710 

0.003 


Table 8.2: Data of Figure 8.8 including low densities. 
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the Fermion sign and the additional spin variable from having two loops instead of one: 



In the spin-polarized case, the exchange diagram entirely cancels the two pair bubbles for 
densities low enough and as a consequence RPA-I-SOSEX excatly cancels in the limit /cp —?• 0. 
RPA-I-APX becomes even positive, since APX contains more diagrams than RPA. This occurs, 
however, only at very low densities beyond ~ 88. 


8.3 APX for solids 

In solids or in molecules the Adjacent Pairs Exchange correction can be evaluated in a similar 
fashion as the Random Phase Approximation. Collecting positive and negative imaginary 
frequencies and rotating the matrices cyclically converts all involved matrices into symmetric, 
real valued matrices; 

vi (l - V^Xo(m)V^) | . 

This makes the evaluation of the inverse numerically more stable. The independent parti¬ 
cle polarizability Xq and the exchange polarizability can be evaluated in imaginary 

frequency according to (6.10) and (8.16) in 0{N^) and in 0{N^) steps, respectively. The 
memory requirement for both cases is 0{N‘^). Accepting a memory footprint of 0{N^), the 
evaluation of the APX and, similarly, of the AC-SOSEX correction can be improved to 0{N‘^) 
by storing the intermediate tensor r(ij^) 

roo 

Pxx'x"(^^) / dv e Gqxx'(^^) ^Ox'x" ( 1^) 

Jo 

and then evaluating the exchange polarizability in terms of E: 

Px^\iX2(i^) = ~ [[ dx3dx4 I - ^ -rrx4XiX3(+ij^)rx3X2X4(-ii^) • 

J J Ts - r4| 


8.3.1 /c-points and Gmax convergence 

In the diagrams of the Random Phase Approximation each Coulomb interaction mediates 
the same momentum q according to momentum conservation at each vertex. This makes the 
RPA the most important contribution for small momenta or, equivalently, at long distances. 
The diagrams of the Adjacent Pairs Exchange correction contain one Coulomb interaction 
with a different momentum than all others. Thus, the APX correction is more short ranged 
than the Random Phase Approximation resulting in a faster /c-point convergence, as shown 
in the inset of Eigure 8.9. 
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Figure 8.9: Basis set extrapolation of the APX energy for LiF with 10 electrons. The black 
points show the energies obtained at different cutoff energies Fimax- The colored points show 
a fit using one, two or three orders in the expansion of the APX energy with respect to large 
cutoff momenta. Each fit only uses data points with a cutoff energy below or equal to that of 
the fitted point. The inset shows the convergence of the APX energy with respect to different 
k point meshes. 

The downside of the short ranged nature of the APX is that the convergence of the 
APX energy is slower with respect to the highest momentum Gmax contained in Xq and 
Assuming that the system is homogeneous at distances short enough, the asymptotic 
behavior of the missing APX energy is equivalent to that of the Uniform Electron Gas, shown 
in Subsection 8.2.1. Eigure 8.9 shows the convergence of the APX energy with respect to the 
employed energy cutoff Emax = Note that the automatic finite basis set extrapolation 

implemented in VASP cannot be used for the APX correction for two reasons. First, unlike 
Xq, Px^^ changes with each considered cutoff momentum since the Coulomb kernel changes. 
Recalculating the exchange polarizability is, however, too time consuming. Second, it is often 
not sufficient to consider only the leading order term of the asymptotic behavior. For an 
accurate extrapolation to Gmax —>■ oo it may be necessary to include the second, or higher 
order terms of the large momentum expansion for the Uniform Electron Gas. The asymptotic 
behavior can then be fit to 

^max ^max V^max 

where ENMAX = G^,^/2 is manually risen. Only the first result line with the largest au¬ 
tomatically chosen cutoff should be used and the number of bands NBANDS should be close 
to the maximum number of bands specified in the OUTCAR file of the DFT or Hartree-Eock 
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Figure 8.10; Relative error of the imaginary frequency integration for different approximations 
of the correlation energy of LiF, depending on the number of samples. The system contains 
10 electrons andjs 384 bands. The inset shows the error for the difference of the APX energy 
for different unit cell volumes, once at 15.4 and once at 15.9 A^. The error is below 1/reV 
beyond the shown number of integration points. 

calculation for each respective value of ENCUT. Note that NBANDS must be a multiple of the 
number of cores employed for the APX calculation. To avoid aliasing effects it is therefore 
advantageous to choose ENMAX such that the reported maximum number of bands is as close 
as possible to NBANDS. 


8.3.2 Frequency grid 

Directly evaluating (8.8) would require 0{N^) steps. However, can be calculated in 

imaginary time using the imaginary time propagators given in (6.11). Subsequently, it can be 
Fourier transformed numerically to imaginary frequency on a non-equidistant grid, analogous 
to calculating the independent particle polarizability as proposed by (Kaltak, Klimes, and 
Kresse 2014b). 

Px^\ix2(i^) = ~ [[ dx3dx4 I- - -1 [ dTie+'''^iGox3Xi(iD)GoxiX4(-iD) 

J J Ts — r4| JO 

poo 

/ dT2e"'''^2Qp^^^^(i.r-2)Gox2X3(-i'^2) (8.16) 

Jo 

The imaginary time propagators can be calculated in 0{N^) steps. If the number of samples 
for the two imaginary time and the final imaginary frequency integration is independent of 
the system size, the evaluation of the APX scales like 0{N^). The prefactor might, however, 
be considerable. 
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In the Random Phase Approximation the employed quadrature frequencies and weights 
are determined from a fit to the function of the direct MP2 energy, since it resembles the RPA 
energy function and its exact frequency integral is known. In the case of the APX energy the 
respective lowest order function would be the exchange MP2 energy 

1 /■ ^ 2(£a - £i){eb - Ej) + 2i/^ _ _1_ 

2 y 27r {{sa - {{eb - SjY + z/2) + Sj - Sa - £b 

and the quadrature frequencies and weights could be fit such that the dominant terms 
with a = b and i = j are best reproduced by the numeric integral 

£j + £^ £a £a ^ ^ ((Fq - EiR + 

for all single particle excitation energies £a — £i- In general, this yields a different grid of 
optimal frequencies and weights than the grid obtained for the Random Phase Approxima¬ 
tion according to (6.12). However, the APX energy also contains the independent particle 
polarizability in all orders beyond the second order, which suggests that the a frequency grid, 
optimized for the exchange MP2 term, might no longer be optimal for higher orders anyway. 
Figure 8.10 shows the convergence of the APX energy in LiF containing 10 electrons with 
respect to the number of imaginary frequency points, using a grid which is optimized for 
the RPA rather than for the APX. The convergence hardly differs for the different energies, 
justifying the use of the RPA optimized grid for calculating and subsequently the APX 

energy. For metallic systemy, knowledge about asymptotic behavior of the APX with respect 
to large imaginary frequencies is important since the quality of the frequency grid degrades 
with vanishing band gap. APX has a different behavior for large frequencies in the Uniform 
Electron Gas than the Random Phase Approximation, as discussed in Subsection 8.2.2. How¬ 
ever, since APX requires a different calculation setup with less /c-points but higher plane wave 
cutoff momenta, evaluating APX and RPA on different frequency grids in metallic systems 
does not pose a considerable disadvantage. 

8.3.3 Lattice constants 

Calculating lattice constants beyond the Random Phase Approximation requires sub-meV 
accuracy for the energy difference at different volumes of the system under consideration. In 
this case, the basis set extrapolation of the APX energy differences does not require more 
than one term in the expansion of the large momentum behavior. In fact, the quality of 
the extrapolation deteriorates when using more than one term in the expansion since the 
higher order terms mostly cancel in the energy difference. Figure 8.11 shows the basis set 
convergence of the energy difference for LiF with 10 electrons. The volumes of the primitive 
cells are 15.9 and 15.4 A^. The data contains considerable aliasing effects since the number 
of bands must be a multiple of the number of cores employed by the calculation, such that 
the number of bands is not always equally close to the maximum number of plane waves 
supported by the respective cutoff energy. All calculations were conducted with a 3 x 3 x 3 
fc-points mesh, providing sufficient accuracy, as indicated in the inset of Figure 8.11. 

At each volume, the Adjacent Pairs Exchange correction is calculated and converged with 
respect to the cutoff energy using VASP with a PBE-DFT reference. The Exact-Exchange-|- 
RPA results were taken from the Birch Murnaghan fit of previous converged RPA calculation. 
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Figure 8.11: Basis set extrapolation of the APX energy difference at different volumes of LiF 
containing 10 electrons. The volumes of the primitive unit cell are 15.9 and 15.4 A^. The 
black points show the energy differences obtained at different cutoff energies Flmax- Aliasing 
effects can be clearly seen. The red points show a fit using in the leading order of the 
expansion of the APX energy with respect to large cutoff momenta. Each fit only uses data 
points with a cutoff energy below or equal to that of the fitted point. Higher order fits offer 
no improvement since they largely cancel in the difference. The inset shows the convergence 
of the APX energy difference with respect to different k point meshes. 
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1000 


Solid 

ao [A] 

RPA 

% 

SOSEX 

% 

APX 

% 

exp. 

C 

(A4) 

3.572 

0.5 

3.552 

<0.1 

3.550 

-0.1 

3.553 

LiH 

(Bl) 

3.983 

0.1 

3.989 

0.3 

3.992 

0.4 

3.979 

LiE 

(Bl) 

3.998 

0.7 

3.955 

-0.4 

3.974 

<0.1 

3.972 

LiCl 

(Bl) 

5.074 

<0.1 



— 

— 

5.070 

NaE 

(Bl) 

4.625 

0.9 



— 

— 

4.582 

MgO 

(Bl) 

4.225 

0.9 



— 

— 

4.189 

GaP 

(B3) 

5.442 

<0.1 



— 

— 

5.439 


Table 8.3: Lattice constants from the Random Phase Approximation, SOSEX and the Adja¬ 
cent Pairs Exchange correction compared to experiment. The experimental lattice constants 
were extrapolated to T = OK and to exclude the zero point energy (ZPE). 
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employing the same PAW potentials. The resulting RPA+APX energy is then fit to a Birch- 
Murnaghan equation of state to retrieve the RPA+APX lattice constants. The pressure 
derivative of the bulk modulus was taken from ... Table 8.3 lists lattice constants from RPA, 
SOSEX and APX in comparison with experiment, extrapolated to T = OK and corrected to 
exclude the phononic zero point energy. 

8.4 Summary and discussion 

The Random Phase Approximation exhibits a systematic error which origins, at least partially, 
from violations of the Pauli exclusion principle. As discussed in the beginning of Chapter 7, 
violating contributions are canceled by diagrams where the propagators of the offending states 
are exchanged. In the lowest order only one pair of propagators is exchanged. 

The Adjacent Pairs Exchange (APX) correction to the RPA constitutes the largest set of 
Eeynman diagrams correcting as many violations as possible in lowest order without intro¬ 
ducing new violations and still forming one ring only: 



This makes the Random Phase Approximation plus the Adjacent Pairs Exchange correction 
a purely ab-initio choice of diagrams when the memory scaling for their evaluation must not 
exceed 0{N‘^), since the ring form of the APX allows its evaluation using only polarization 
parts with 2 open ends and RPA is the most important class of diagrams for high densities, 
becoming exact in the limit rg —)• 0. This argument relies on the prevalence of long ranged 
states making the RPA an accurate approximation in the first place, which might not hold 
in the presence of more localized states at the Permi edge. However, in cases where the RPA 
offers a good approximation, no test against experiment is required to argue that APX offers 
the lowest order correction - which makes it ab-initio'^. 

Table 8.4 compares the SOSEX, the AC-SOSEX and the APX correction to the Random 
Phase Approximation in terms of contained Goldstone diagrams, along with their respective 
computational costs in time and memory. The diagrams of SOSEX constitute the largest set 
of Goldstone, rather than Eeynman, diagrams correcting as many violations as possible under 
the same conditions than in the APX. In the case of SOSEX, forming a ring is necessary for the 
applicability of the direct ring Coupled Cluster Doubles amplitudes, which can be constructed 
in 0{N^) steps, as discussed in Section 6.2. Dropping the condition for the SOSEX diagrams 
to form a ring raises the time complexity for the construction of the amplitudes to 0{N^). 
Accepting this time complexity and the memory complexity of SOSEX, the full Coupled 
Cluster Singles and Doubles (CCSD) approximation can be employed. It includes considerably 
more diagrams, as indicated in Table 8.4, relieving the need of a correction to the Random 
Phase Approximation in general. In cases, where CCSD is known to be less accurate, as 
for instance in dissociation processes, the distinguishable Coupled Cluster method may offer 
improvement at the same computational costs. 

^Still, only a test against the uniform electron gas induced further examination of the two sided APX and 
the possibility of new violations introduced by it. Testing against systems - theoretical or not - is therefore 
necessary, of course. 
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Table 8.4; Comparison of different approximations beyond the Random Phase Approximation, 
showing the lowest order Goldstone diagrams introduced by the respective approximation. 
The AC-SOSEX is not derived within the same many-body perturbation theory framework 
as the other approximations. It can, however, be translated into Goldstone diagrams when 
including swapped ladder diagrams, as shown here in third order and which are discussed in 
Section 7.3. 

Including exchange diagrams to cancel violating contributions in the RPA diagrams offers 
an ab-initio strategy for selecting the lowest order diagrams of many-body perturbation the¬ 
ory, correcting for RPA’s systematic error. In the framework of the Adiabatic Connection, the 
unknown exchange-correlation kernel occurring in the Dyson-like equation (6.27) of the 
polarizability X;,, has to be included for a correction to the RPA. The AC-SOSEX contains 
the exchange-correlation kernel in lowest order: 

vAC—SOSEX -v" fA "Y" I "v pA -v" \A7”Y" i 

= AoIxc-^0 + v^Olxc^OA V + • • • , 
approximating the kernel by which is given implicitly by 



In the first and in the last diagram, particles turn into holes and vice versa at the Coulomb 
interaction, which has no correspondence in many-body perturbation theory. This is, there¬ 
fore, not an obvious choice for an approximation in the author’s opinion. It neither originates 
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from time dependent DFT nor from many-body perturbation theory, which is employed by 
AC-RPA after all to calculate the independent particle polarizability Xq. As discussed in Sec¬ 
tion 7.2, the above approximation rather comes from transforming the RPA energy expression 
into the particle/hole basis, applying the analogy of SOSEX to exchange two indices, ignoring 
the time order, and then transforming it back to the position basis. However, the results are 
very similar to those of SOSEX for reasons discussed in Section 7.3, which are in no obvious 
connection to its derivation. Owing to its similar results to SOSEX, the AC-SOSEX can be 
considered a numerical recipe to approximate the SOSEX energy, requiring only 0{N‘^) of 
memory. 

The Adjacent Pairs Exchange corrections contains more diagrams than the Second Order 
Screened Exchange correction since APX is the pendant of SOSEX in terms of Eeynman 
diagrams rather than Goldstone diagrams. This does, however, not imply that APX is more 
accurate in any given situation. After all, the surprisingly high accuracy of SOSEX or APX 
in case of the uniform electron gas is certainly fortuitously and no theory, involving only first 
and second order exchange diagrams can be expected to provide this accuracy in general. 
The spin-polarized uniform electron gas is an example of a system, where RPA-I-SOSEX or 
RPA-|-APX can yield zero or even positive correlation energies in the limit of dilute densities. 
Por densities of real metals, the performance in the spin-polarized case is better. APX misses 
an accurate correction to RPA by less than 20% and SOSEX is slightly better. It is assumed 
that AC-SOSEX is close to SOSEX since no SOSEX calculations for the spin-polarized UEG 
in the thermodynamic limit have been found. It is worth mentioning that finite spin-polarized 
systems often exhibit a (quasi)-degenerate ground state, rendering many-body perturbation 
theory, as discussed here and employed by most implementations, inaccurate at best. (Quasi)- 
degenerate many-body perturbation theory is discussed, for instance, in (Shavitt and Bartlett 
2009). 

APX and SOSEX are identical up to third order resulting only in small differences. In the 
case of the non-spin-polarized uniform electron gas, RPA-I-SOSEX matches Quamtum Monte 
Carlo correlation energies already at Vg ~ 5, where RPA-I-APX still has an error of about 
than 0.7 mE^ per electron, which is a relative error 5%. Eor lower densities with Vg > 8, 
where correlation effects are more prevalent, APX improves on SOSEX and its accuracy is 
never worse than 0.3mE/j per electron, even up to Vg = 50. Note, that the accuracy is given 
with respect to the (Perdew and Zunger 1981) fit of the (Ceperley and Alder 1980) QMC 
resutls. The Adjacent Pairs Exchange correction can be computed in 0{N^) with a memory 
requirement scaling like 0{N‘^). This is equivalent to the AC-SOSEX but considerably less 
memory demanding than SOSEX. If a memory demand of 0{N^) is permissible, APX and AC- 
SOSEX can also be evaluated in 0{N‘^) steps. Both, APX and AC-SOSEX are implemented 
in VASP and first applications on lattice constants show excellent agreement with experiment. 
A more extensive survey including atomization energies is still to be made. 

Not all improvements beyond the Random Phase Approximation are based on the inclusion 
of exchange processes. In the GW approximation the iterative scheme of RPA, used to build 
an effective interaction W from the Coulomb interaction V according to (6.7), is extended 
by an iterative scheme to improve the description of the propagator Go towards an effective 
propagator G of the fully interacting system, referred to as dressed propagator. There is a 
whole family of such approximations, based on the work of (Hedin 1965). The central quantity 
of interest in these approximations is the dressed propagator, rather than connected diagrams, 
as in the Goldstone approach of many-body perturbation theory, discussed here. This makes 
a comparison between the two approaches hard, as there is no consideration of symmetries 
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and the propagators are approximated by merely shifting the poles of the initial propagators 
Gq. In most cases the dressed propagator is used to retrieve excitation spectra but there are 
also total energy calculations within the GW approximation, usually employing the formula of 
(Galitskii and Migdal 1958). In Appendix A, an alternative approach is suggested to evaluate 
the total energy in the GqWq approximation strictly within the many-body perturbation 
theory discussed here. The total energies retrieved by this approximation are expected to 
be more accurate than those of the RPA at computational costs that should not exceed 
0{N‘^) but this remains to be tested in the future. Also, GqWq is only the least accurate 
approximation of the mentioned family of approximations, depending on the reference such 
as DFT or HF and further, exchange effects are known to be important not only to correct 
for violations of the exclusion principle, so it is unclear, whether GqWq total energies are a 
viable option to RPA+APX or RPA+AC-SOSEX. 

Finally, it is worth mentioning that most of time and memory requirement of high accu¬ 
racy methods, such as Coupled Cluster Singles and Doubles (CCSD) strongly depends on the 
number of unoccupied states ipa, needed for a convergent result. This number can be reduced 
without considerably sacrificing the accuracy such that CCSD or even higher accuracy meth¬ 
ods become feasible for larger systems. One way of reducing the number of unoccupied states 
is by means of natural orbitals, applicable in the case where large voids are between the atoms 
or molecules, such as in atomization energy calculations of solids. Another way of reducing 
the number of unoccupied states is by including explicit correction already in the descrip¬ 
tion of the unperturbed system. This can be done by augmenting the Slater determinant, 
which is a product of functions depending on one electron position only, by a set of functions 
/(xi,X 2 ) explicitly depending on two electron positions, hence the name /12 methods. In all 
cases there is still an underlying perturbation expansion to be evaluated and the Adjacent 
Paris Correction as well as AC-SOSEX can also profit from such a reduction of the number 
of unoccupied states. 
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Appendix A 

Total energies in GqWq from 
connected diagrams 


Given the propagator Go in matrix form and in imaginary time according to (6.11)^ 

{ - forr<0 

i 

+ ^ otherwise, 


Goxx'(i'r) = { 


we can Fourier transform it to imaginary frequency Go(ir/)- We also write the screened 
Coulomb interaction Wq in the Random Phase Approximation (RPA) in the same form 


w 

— 


+ 


Wo(R) = 

V 

+ 




+ 




+ 


(A.l) 


= V 1 - XnfmlV 


-1 


where Xq is the independent particle polarizability in imaginary frequency as retrieved from 
a Fourier transform of (6.10). The irreducible self energy S can then be approximated in the 
RPA by 


:a*0 


+ 


(A.2) 


idz/ 


51oxx'(W) = / i^Goxx'(i^ + i«^) Woxx'(iz^) - Vxx' • 


Note that the matrices are multiplied elementwise and that we exclude the exchange term GqV 
since it cancels with the effective interaction in a Hartree-Fock reference as discussed in Section 

^Note that the propagator used here differs by a factor of i from the usual definition of the propagator, as 
discussed in the footnote of (5.42) 
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APPENDIX A. TOTAL ENERGIES IN GqWo FROM OONNECTED DIAGRAMS 


5.6. In the GqWq approximation the full propagator G is approximated by inserting the above 
approximation to the irreducible self energy Sq into the propagator of the unperturbed system 
Go arbitrarily many times 


G 


GnSnGr 


Go^IoGo^IoGc 


where we omit the imaginary frequency argument for brevity. Properties of the interacting 
system can then be extracted from this approximation to the full propagator G. To get the 
total energy the formula of (Galitskii and Migdal 1958) can be used on the full imaginary 
time propagator G(ir), retrieved from an inverse Fourier transform of G(i? 7 ) at r —)■ 0 from 
above: 

E = -^ [ dx lim lim ( Grar'ai^r), (A.4) 


dx lim lim 

r'—>-r r—>-0+ 


where /ix is the single body Hamiltonian of the unperturbed system acting on Gx. This 
can be interpreted as cutting out one occurrence of Go by the differential operator and then 
closing the remaining diagrams contained in G (Ziesche 2010). 

Here, an alternative way to the Galitskii-Migdal formula is proposed for evaluating the 
total energy in GoVFo) respecting the symmetries of all connected diagrams occurring in this 
approximation, as discussed in Subsection 5.8.2. Figure A.l shows all diagrams order by 
order, where each row contains an order of the expansion of Wq contained in Sq according to 
(A.l) and each column contains an order of G according to (A.3). Note that the symmetries 
of the first column differ from the symmetries of all other columns since closing Sq with one 
Go forms a bubble equivalent to Xq. The diagrams of the first column form the ring diagrams 
of the RPA. All the remaining diagrams have rotational symmetry but no reflection symmetry 
and can be summed to infinite order in the same fashion as it was done in the RPA in (6.8), 
arriving at 


^GoM/o^^RPA^i I ^t,|_l('G^(i^)So(ir?))'-^(Go(ir/)So(ir?))'-...| 

= Ef^^ - j ^ tr| log (l - Go(i??)So(i7?)) + Go(ir/)Xo(i? 7 )} . (A.5) 

This approach still has to be tested and compared to previous GqIPo total energies, as for 
instance retrieved by (Holm and Aryasetiawan 2000). The two approaches do differ after all. 
Gonventional Go Ho uses a normalization factor for calculating G while this approach uses 
symmetry factors. An implementation for the uniform electron gas in the thermodynamic 
limit is still ongoing. 
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Figure A.l: Connected (closed) diagrams occurring in the GqWo approximation of the cor¬ 
relation energy. The parts of the self energy are given in blue while the parts of the closing 
propagator are given in red. With each column the number of insertions of SqGo into G 
increments, while the number of insertions of XqV into So increments with each row. The 
single Fermion loop in the center results in a negative Fermion sign of all diagrams beyond 
RPA. Imaginary units, the imaginary frequency integration and the trace are omitted for 
clarity. 









